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Zusammenfassung 



Diese Thesis umfasst mehrere Aspekte theoretischer Stellardynamik in Sternhaufen, sowohl 
in analytischer als auch in numerischer Hinsicht. Wir versuchen, Licht auf Phanomene zu 
werfen, welche zur Zeit in alien Galaxietypen beobachtet werden, einschlieBlich AGNs und 
Quasare, welche zu den machtigsten Objekten des Universums zahlen. Die Wechselwirkun- 
gen zwischen einem Sternsystem und einem zentralen schwarzen Loch fuhren zu einer Menge 
interessanter Phanomene. Die von uns verwendete Methode ermoglicht eine Betrachtung 
leicht einsichtiger Aspekte ohne jegliches Rauschen, welches die Teilchen-Methoden mit 
sich bringen. Wir untersuchen die wichtigsten physikalischen Prozesse, die in der Entwick- 
lung eines spharischen Sternhaufens ablaufen, etwa Selbstanziehungskraft, Zwei-Korper- 
Relaxation etc sowie die Wechselwirkung mit einem schwarzen Loch und die Funktion des 
Massenspektrums. Wir beschaftigen uns jedoch nicht nur mit diesem Thema alleine, sondern 
auch mit einer Analyse supermassiver Sterne. Wie diese Sterne die Aktivitaten der Quasare 
durch Sternakkretion und Energiestrom antreiben ist eine der Fragen, die hierbei aufkom- 
men. Wir gehen auch anderen Fragen nach, etwa jener nach der noch nicht verstandenen 
Entwicklung eines solchen Objektes und seiner Wechselwirkung mit dem ihn umgebendem 
Sternsystem. Dies ist ein Kernpunkt der Astrophysik, da diese Objekte als Vorlaufer von 
supermassiven schwarzen Lochern betrachtet werden konnen. 



Abstract 

This thesis embraces several aspects of theoretical stellar dynamics in clusters, both analyti- 
cally and numerically. We try to elucidate the phenomena currently observed in all types of 
galaxies, including AGNs and quasars, some of the most powerful objects in the universe. 
The interactions between the stellar system and the central black hole give rise to a lot of 
interesting phenomena. The scheme we employ enables a study of clean-cut aspects without 
any noise that particle methods suffer from. We study the most important physical processes 
that are readily available in the evolution of a spherical cluster, like self-gravity, two-body 
relaxation etc, the interaction with a central black hole and the role of a mass spectrum. Not 
only embark we upon this subject, but we set about an analysis on super-massive stars. How 
these stars could power the quasar activity by star accretion and energy flows is one of the 
questions that arises. We undertake other questions, such as the uncertain evolution of such 
an object and its interaction with the surrounding stellar system. This is of crucial importance 
in astrophysics, for these objects could be regarded as super-massive black holes progenitors. 



Antonio Amaro Pita, my father, never could see finished this work, because he passed away in Novem- 
ber of last year. It is very difficult to write these lines, because I miss him a lot. He lived intensively, 
much more than all people I know. He squeezed the nice things out of life until the very last drop and 
then he decided to go away and carry on with the fun somewhere else because here it was damned 
boring for him. 

When my father died, Marc Freitag, who was in the north of America, wrote: 

En Pennsylvanie c 'est un automne magnifique avec des arbres qui deviennent tres rouges avant de 
se depouiller totalement. J 'imagine qu'il est normal que les grands arbres, devenus vieux, soient 
terasses par le vent... Cela n'enleve rien a la beaute de laforet. 

There is death because there is life, and there is life because there is death. My father died, and so 
will I someday; this is not a reason to be upset, it is just the flowing of life. It steals us the surface of 
the sand and allows us to walk in the new surface and leave our prints. Sabine is now pregnant and I 
hope our child will be as full of life as his grandfather was and able to enjoy every second of his new 
life, just as my father did. 
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hen he was in this plight, Ino daughter of Cadmus, also called Leu- 
cothea, saw him. She had formerly been a mere mortal, but had been 
since raised to the rank of a marine goddess. Seeing in what great 
distress Ulysses now was, she had compassion upon him, and, rising 
like a sea-gull from the waves, took her seat upon the raft. "My poor 
good man," said she, "why is Poseidon so furiously angry with you? 
He is giving you a great deal of trouble, but for all his bluster he will not kill you. You 
seem to be a sensible person, do then as I bid you; strip, leave your raft to drive before 
the wind, and swim to the Phaecian coast where better luck awaits you. And here, take 
my veil and put it round your chest; it is enchanted, and you can come to no harm so 
long as you wear it. As soon as you touch land take it off, throw it back as far as you 
can into the sea, and then go away again." With these words she took off her veil and 
gave it him. Then she dived down again like a sea-gull and vanished beneath the dark 
blue waters. But Odysseus did not know what to think. "Alas," he said to himself in his 
dismay, "this is only some one or other of the gods who is luring me to ruin by advising 
me to will quit my raft. At any rate I will not do so at present, for the land where she 
said I should be quit of all troubles seemed to be still a good way off. I know what I will 
do- 1 am sure it will be best- no matter what happens I will stick to the raft as long as 
her timbers hold together, but when the sea breaks her up I will swim for it; I do not see 
how I can do any better than this. 

Odyssey, Homer (Book V) 




Odysseus clearly distrusted the benevolence of the goddess. And his resolution was quite clear: neither 
to accept the advice nor the favours of superhuman origin until tragedy be imminent. While the ship re- 
mained whole, he would not dare leave it... despite the clear help that (whichever) gods offered. Perhaps 
because Leucothea had once been mortal? This did not seem to be the reason, in my view. Odysseus, 
one of the first western heroes of whom we have written evidence, did not rely at all on the gods who 
favoured him. 

Alone, wet, cold, in a feeble boat that could not guarantee him the necessary cover against the wrath of 
another god, in the middle of the night, battling the fearsome waves that carried him directly to the reefs, 
Odysseus looks at the piece of sail that the goddess gave him with which he will have to reach the cost, 
still far away. Perhaps because he is a hero, he arrives there where many do not. Those who, desperate, 
would blindly have heeded the advice of gods. He, on the contrary, proceeds rationally. Odysseus will 
not pay attention to the advice of gods until his chances be clearly doomed, until his ship disintegrates, 
taken by the rage of the sea, and he has nothing to lose. Then, and only then, he will abandon the chunk 
of wood in which he finds himself at the mercy of the storm, and he would launch himself into the water 
with the hope that the madness of the seagull-goddess be true. 

Naturally, as the reader can guess, either because he is familiar with Homer's work or sufficiently savvy 
to be irrational and see it, Odysseus reaches the shore. 

I read this passage a couple of years ago, when I was still studying in Valencia. It impressed me to see 
that a mythological greek hero, maybe the mythological greek hero could think in such a rational way 
and only leave his last shelter when all other possibilities were ruled out. I interpreted it like a knowing 
wink of Homer to us, poor mortals. When I started my PhD I felt a bit like Odysseus. Only a couple of 
months ago I decided to trust the veil and now, safe in the cost, I resolved to include this text here. 
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Motivation 



1.1 What is this all about? First words 

MASSIVE objects avoiding light to escape from it is a concept that goes back to the 18th century, 
when John Michell (1724 - 1793), an English natural philosopher and geologist overtook 
Laplace with the idea that a very massive object could be able to stop light escaping from it 
thanks to its overwhelming gravity. Such an object would be black, that is, invisible, precisely because 
of the lack of light (Michell, 1784; Schaffer, 1979). He wrote: 

"If the semi-diameter of a sphere of the same density as the sun is in the proportion of five 
hundred to one, and by supposing light to be attracted by the same force in proportion to its 
mass with other bodies, all light emitted from such a body would be made to return towards 
it, by its own proper gravity." 

Even thought we shall sift through this concept oftentimes in the work to be presented along the next 
chapters, it appears to be necessary to give an initial, short overview of the problem's subject we shall 
tackle in this thesis. 

A "black hole' embraces the observation of phenomena which are associated with matter accretion on 
to it, for we are not able to directly observe it. Emission of electromagnetic radiation, accretion discs 
and emerging jets are some, among others, evidences we have for the existence of such massive dark 
objects, lurking at the centre of galaxies. 

On the other hand, spectroscopic and photometric studies of the stellar and gas dynamics in the inner 
regions of local spheroidal galaxies and prominent bulges give us the certainty that nearly all galaxies 
should harbour a central massive dark object, with a tight relationship between its mass and the mass or 

'This term was first employed by John Archibald Wheeler (b. 1911), an American theoretical physicist 
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the velocity dispersion of the host galaxy spheroidal component (as we will see ahead). Nonetheless, 
even though we should mention that we do not dispose of any direct evidence that such massive dark 
objects are black holes, alternative explanations are sorely constrained (see, for instance, Kormendy 
2003). 

Super-massive black holes are ensconced at the centre of active galaxies. What we understand with ac- 
tive is a galaxy in which we can find an important amount of emitted energy which cannot be attributed 
to its "normal" components. AGNs, as they are usually denominated, have the powerhouse for their 
observed phenomena in a compact region in the centre. 

We will embark in the next chapters on an analytical and numerical study of the dynamics of stellar 
systems harbouring a central massive object in order to extract the dominant physical processes and 
their parameter dependences like dynamical friction and mass accretion. 

This chapter is devoted to a general description of the scenario with which we shall deal with. Firstly, 
we give a short bird's-eye view of active galactic nuclei in order to clearly exhibit the boundaries of 
our problem and then, secondly, we will give an extending to the present time summary on the possible 
nature of the central dark object and its possible origin and formation. 

1.2 What are AGNs and what makes them interesting? 

The expression "active galactic nucleus" of a galaxy (AGN henceforth) is referred to the energetic phe- 
nomena occurring at the central regions of galaxies which cannot be explained in terms of stars, dust 
or interstellar gas. The released energy is emitted across most of the electromagnetic spectrum, UV, 
X-rays, as infrared, radio waves and gamma rays. Such objects have big luminosities (10 4 times that 
of a typical galaxy) coming from tiny volumes (<C lpc 3 ); in the case of a typical Seyfert galaxy the 
luminosity is about ~ 10 11 L & (L© = 3.83 • 10 33 erg/s is the luminosity of the sun), whilst for a typical 
quasar it is brighter by a factor 100 or even more; actually they can emit as much as some thousand 
galaxies like our Milky- Way. They are thus the most powerful objects in the universe. There is a con- 
nexion between young galaxies and the creation of active nuclei, because this luminosity can strongly 
vary with the red-shift. 

Anticipating something on which we will elaborate a bit ahead, nowadays one explains (although this 
is still not accepted ubiquitously 2 ) the generation of energy resorting to matter accreting on to a super- 
massive black hole in the range of mass <~ 1O 6 ~ 1O M (where is the black hole mass). In this 
process, angular momentum flattens the structure of the in-falling material to a so-called accretion disc. 
It is frequent to observe jets, which may be arising from the accretion disc (see Fig. 1.1), although we 
do not dispose of direct direct observations that corroborate this. Accretion is a very efficient channel 

2 For some alternative and interesting schemes, see Ginzburg and Ozernoy (1964) for spinars, Arons et al. (1975) for clusters 
of stellar mass BHs or neutron stars and Terlevich (1989) for warmers: massive stars with strong mass-loss spend a significant 
amount of their He-burning phase to the left of the ZAMS on the HR diagram. The ionisation spectrum of a young cluster of 
massive stars will be strongly influenced by extremely hot and luminous stars 
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of turning matter into energy. Whereas nuclear fusion reaches only a few percent, accretion can breeze 
in almost 50% of the mass-energy of a star into energy. 

Being a bit more punctilious, we should say that hallmarks for AGNs are the frequency of their electro- 
magnetic emission frequencies, observed from < 100 MHz (as low frequency radio sources) to > 100 
MeV (which corresponds to <~ 2 • 10 22 Hz gamma ray sources). Giant jets give the upper size of manifest 
activity < 6 Mpc ~ 2 • 10 25 cm 3 , and the lower limit is given by the covered distance by light in the 
shortest X-ray variability times, which is <~ 2 • 10 12 cm. 

As regards the size, and resorting again to Fig. (1.1), we can envisage this as a radial distance from the 
very centre of the AGN where, ostensibly, a SMBH is harboured, and the different observed features of 
the nucleus. From the centre outwards, we have first a UV ionising source amidst the optical continuum 
region. This, in turn, is enclosed by the emission line clouds (NLR, BLR 4 ) and the compact radio 
sources and these betwixt another emitting region etc. 

The radiated power at a certain frequency per dex 5 frequency ranges from ~ 10 39 erg/s (radio power 
of the MW) to <~ 10 48 erg/s, the emitted UV power of the most powerful, high-redshifted quasars. 
Such broad frequencies and radius ranges for emission makes us to duly note that they are far out of 
thermal equilibrium. This manifests in two ways: first, smaller regions are hotter; second, components 
of utterly different temperature can exist together, even though components differ in one or two orders 
of magnitude in size. 

1.2.1 Outstanding features of AGNs 

We have to somehow specify what we mean with the term AGN. For this scope we should name the 
observable phenomena which are used to detect them; on the other hand, nevertheless, AGNs can be 
noticed in many ways, but this does not mean that all of them partake the same features. In fact, an AGN 
"picks up" some of these defining qualities to be such an object. 

Very small angular size If we have optical images of the host galaxy, the nucleus happens to be a 
very bright point with a flux that can touch or even surpass the rest one of the host galaxy. In Fig. (1.2) 
we can see an example of this: NGC 1566, a nearby AGN (Jarrett et al., 2003). 

Radio astronomical observations were the first proof for the existence of non-stellar activity in ex- 
ternal galaxies. Early observations revealed that many bright radio sources have the shape of two lobes 
with a galaxy situated at half the distance (see Fig. 1.3). Even though radio astronomical techniques 
are powerful, they might be misguiding, since the radio band never accounts for more than <~ 1% of 

3 If we do not take into account the ionising radiation on intergalactic medium 

4 See subsection 1.2.2 

5 A "unit" in the logarithmic axis 
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14 15 16 17 18 19 

log ]o r (cm) 

Fig. 1.1: Unified model of an AGN. A relativistic jet is only to be found in radio-loud sources. In the vicinity 
of the central BH we may find a flared-up disc due to X-rays. 

the bolometric luminosity; to boot, less biased surveys bore witness that the most AGNs emit a much 
smaller portion of their total power in the radio. 

High luminosity This can run in the case of an AGN from about 1% of a typical galaxy up to ~ 10 4 
times as great. Of course these limits are not completely fixed, for there can be "hidden" a huge popula- 
tion with lower, non observable luminosities; another possible reason for the incapacity of observation 
is relativistic beaming, or obscuration, due to thick dust extinction, which can misguide our measure- 
ments. Compared with the galaxy spectra, an AGN continuum spectra looks stunningly different. The 
energy flux per logarithmic bandwidth produces a spectrum far broader compared with that of a galaxy. 
In comparison with a galaxy, AGNs give out in the radio band a fraction of bolometric luminosity that is 
about an order of magnitude greater, and the corresponding in the x-ray band is between three and four 
larger. 
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Fig. 1.2: NGC 1566 (Jarrett et al., 2003) 




Fig. 1.3: Cygnus A at 6 cm wavelength (Perley et al., 1984). The total extent of this source is ~ 120 kpc 



Variability is often said to be a distinctive characteristic of AGNs; this is relative for, even though 
the most of them can be seen to vary in the optical band, the typical amplitude over time-scales is 
frequently 10%. There is no a fixed time-scale on which AGNs vary, and therefore it is hard to measure 
their amplitude of variability. We can find a subgroup of these objects that can be observed to vary even 
from night to night, and cumulative changes of factors of hundred have happened over year time-scales. 
In Fig. (1.4) we can see an example of this rapid variation (Pian et al., 1997). 
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20 22 
May 1994 



Fig. 1.4: Two ultraviolet light curves for PKS 2155-304. The open circles show the continuum flux at 2800 
A; the filled ones at 1400 A. We can see that within one day changes of several tens of percent occurred 



Polarisation Like in the case of galaxies, even though stars are in itself unpolarised, the light we 
observe is by and large polarised by about 0.5% because of interstellar dust transmission polarisation. 
As for the AGNs, they are also polarised, linearly and with a fractional polarisation of between 0.5- 
— 2%. We find a minority that is much more polarised, frequently ~ 10%. We have a strong variation 
on both the magnitude and the direction of the polarisation for those objects that are strongly polarised 
and have strongly variable total flux. Nonetheless, variability does have bounds, for circular polarisation 
has never been detected. 



Emission lines This draws our attention, for it is easy and productive to study them. On the contrary 
to most of galaxies, these are very prominent for this kind of objects. The observed emission lines are 
stereotypic from one AGN to the other; almost at all times we observe Lya, the Balmer lines, the Civ 
1549 doublet etc. Oftentimes is seen the FeKa x-ray line near 6.4 KeV. 
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1.2.2 AGN taxonomy 

To key out the sub-variety of AGNs we are talking about we have several terms. Whether or not the 
distinctions is beneficial is something we could call into question, but we have to know it if we want 
to follow "the conversation". A deep description of all AGNs nomenclatures with their typical charac- 
teristics is beyond the scope of this work. Here we just want to give a general idea of the vast variety 
of AGNs and main features. We have different ways of categorisation, depending on by which crite- 
ria we want to classify. Some of the terms are self-explanatory; like, for instance radio loud or radio 
quiet. OVV is the acronym for Optically Violently Variable, for in the optic band we have very rapid 
and large amplitude variability, as the name itself reveals. Others are coined terms; Quasar is simply 
the pronounced form of "QSR" (Quasi-Stellar Radio source), but after a time the meaning changed and 
evolved to "generic AGN"; nowadays it has nothing to do with the radio luminosity. In fact, for instance, 
very low luminosity AGNs are called micro-quasars. Some terms bring up the name of the first people 
that identified the class: Carl Seyfert remarked the first Seyfert galaxies, and they were split up later into 
two types, depending on the existence or not of broad wings in the emission lines. In an equal manner, 
Fanaroff and Riley pointed out a distinction in luminosity and morphology among the radio galaxies, 
and so they were named FR galaxies after them. We find also far-out histories in the nomenclature of 
AGNs: The type BL Lac Objects were initially identified as variable stars in the Lacerta constellation 
and, thus, they were named "BL Lac". The term Blazar springs up from the fact that the power output 
of the classes OVV and BL Lac "blazes" dramatically and was thought to unify them, because they are 
very similar. 

In table 1.1 (taken from Krolik 1999) an object is said to be point-like if an optical point source can be 
seen. By broad-band we mean that there is a comparable luminosity in the infrared, optical and x-ray 
bands. The existence in the optical and ultraviolet spectra of lines several thousand (hundred) km/s in 
width is meant with broad-lines (narrow-lines). By radio we understand that for the object the fraction 
of luminosity emitted in the radio is relatively large, perhaps <~ 10~ 3 of the bolometric. The members of 
a class should vary by an order of magnitude or more in the optical band over a human time life for the 
variable entrance and for polarisation the optical light should be at least a few percent linearly polarised. 
In the table we have added the subdivision LINERs, which stands for "Low-Ionisation Nuclear Emission 
Region" (Heckman, 1980), but this category is on the verge of activity; it is not clear whether or not 
these galaxies are really active. 

This classification allows us to arrange the types in a three-dimensional parameter space, since we find 
that we have the groups radio-loud against radio-quiet, strong variable against all the others and narrow 
emission lines against broad emission lines. We can display a three-parameter space in which we locate 
each of the types of AGNs (after Krolik 1999), as we can see in Fig. (1.5). 

1.2.3 The unified model 

After a study of almost 20 years, as we commented at the beginning of this chapter, a unification scheme 
for AGNs has appeared in the community. According to it, a well-mannered AGN should have the 
following identifying characteristics (see Fig 1.1), 
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Type 


Point-like 


Broad-band 


Broad lines 


Narrow lines 


Radio 


Variable 


Polarised 


Radio-loud quasars 


yes 


yes 


yes 


yes 


yes 


some 


some 


Radio-quiet quasars 


yes 


yes 


yes 


yes 


weak 


weak 


weak 


Broad line radio galaxies 


yes 


yes 


yes 


yes 


yes 


weak 


weak 


(FR2 only) 
















Narrow line radio galaxies 


no 


no 


no 


yes 


yes 


no 


no 


(FR1 and FR2) 
















OVV quasars 


yes 


yes 


yes 


yes 


yes 


yes 


yes 


BL Lac objects 


yes 


yes 


no 


no 


yes 


yes 


yes 


Seyferts type 1 


yes 


yes 


yes 


yes 


weak 


some 


weak 


Seyferts type 2 


no 


yes 


no 


yes 


weak 


no 


some 


LINERs 


no 


no 


no 


yes 


no 


no 


no 



Table 1.1: The AGN taxonomy 
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radio 



NLRGs 



BL Lacs 



BLRGs 
RLQs 



Seyfert 2s 



OVVs 



variability 



Seyfpft Is 
RQ^s 
emission line width 
Fig. 1.5: Classification of AGNs in a three-dimensional parameter space 



Super-massive black hole -in the range of 10 6 10 M 

Accretion disc and corona - heated by magnetic and/ or viscous processes with optical radiation though 

soft X-ray energy 
Broad line region - with high velocity gas 
Narrow line region - corresponding to lower velocity gas 

Obscuring medium - which can adopt a torus form or another geometry. This has the characteristic, 
as its names pinpoints, of hiding the broad-line region from some directions owing to its material 
(dust). It is located at about 10 - 100 pc between the smaller inner region from which the broad 
emission lines come from and the more external zone where the narrow emission lines emerge 

Relativistic jet - emerging at distances of about ^100 ^"schw (Schwarzschild radius) of the BH. The 
extension of these objects can be out to tens of hundreds of parsecs or even Mpc 

Excluding intrinsic variances in BH's differentia, like ionisation parameter, size, density, luminosity 
etc, the unification model properties (except for, maybe, relativistic jets) stay staunch to all AGNs. In 
conformity with it, many of the main observational characteristics that we have described before are 
due to orientation and are not intrinsic differences. The galaxy appears as a Seyfert 2 or NLRG at the 
regions where the obscuring medium stops the direct view of the central parts. Facing the AGN, as we 
look into the inner regions and it appears as Seyfert Is or, in radio-loud, BLRGs to radio-loud quasars 
and at inferior luminosities, BL Lacs. 
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1.2.4 Appraising Jt* 

A direct consequence of the paradigm of SMBHs at the centre of ancient galaxies to explain the energy 
emitted by quasars is that relic SMBHs should inhabit at least a fraction of present-day galaxies (Rees, 
1990). This conclusion was first made quantitative by Soltan (1982) and has recently be revisited in 
more detail and in the light of recent observations by Yu and Tremaine (2002). 

In a system of luminosity L with a lifetime 7] um , if e is the conversion of matter to energy (accretion) 
efficiency, one should expect that in the system, an amount of material L7] um /ec 2 has been gathered 
at some place, where c is the light velocity. We are referring to the central super-massive black hole, 
harboured by the host galaxy. We will give a description of such an object in the next section and pore 
over it in the next chapters. 

If we are capable of determining 7] um and e, we will have a rough idea of what M is. Usually, for e one 
can prove that e <~ 1/10 is a suitable value (McCray, 1979). Assuming that the lobes are powered by 
the central engine, we have that, for Cyg A, with a lobe separation of 80 kpc, 7i um <~ 4 • 10 7 yr and, so, 
~ 1O 8 M . If we do not have at our disposal double radio sources, we can reckon in a different way: 
if the lifetime of the galaxy is ~ 10 10 yr and 1% of galaxies are active and they all go through activity 
phases, 7i um ~ 10 8 yr; this means that, for the more luminous sources (L > 10 47 erg/s), ~ 1O 1O M . 
On the other hand, if we had resorted to the luminosity per unit volume from quasars from observations, 
we would have found out that, since this gives the energy output per unit volume over the age of the 
universe, the mass density of "recycled" material is consistent with an average dead mass of 1O 8 M . A 
number of independent deductions leads to a similar result (Frank et al., 2002). 

1.3 Massive black holes and their possible progenitors 
1.3.1 (Super-) massive black holes 

The quest for the source of the luminosities of L rts 1O 12 L produced on such small scales, jets and other 
properties of quasars and other types of active galactic nuclei led in the 60's and 70's to a thorough 
research that hint to the inkling of "super-massive central objects" harboured at their centres. These 
were suggested to be the main source of such characteristics (Lynden-Bell, 1967; Lynden-Bell and 
Rees, 1971a; Hills, 1975). Lynden-Bell (1969) showed that the release of gravitational binding energy 
by stellar accretion on to a SMBH could be the primary powerhouse of an AGN (Lynden-Bell, 1969). 
In the last decade, observational evidences have been accumulating that strongly suggest that massive 
BHs are indeed present at the centre of most galaxies, with a significant spheroidal component. Mostly 
thanks to the HST, the kinematics in present-day universe of gas or stars has been measured in the central 
parts of tens of nearby galaxies. In almost all cases 6 , proper modelling of the measured motions requires 
the presence of a central compact dark object with a mass of a few 10 6 to 1O 9 M (Ferrarese et al., 2001; 
Gebhardt et al., 2002; Pinkney et al., 2003; Kormendy, 2003, and references therein). Note, however, 
that the conclusion that such an object is indeed a BH rather than a cluster of smaller dark objects (like 

6 With, notably, the possible exception of M33 (Gebhardt et al., 2001; Merritt et al., 2001) 
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neutron stars, brown dwarves etc) has only been reached for two galaxies. The first one is the Milky 
Way itself at the centre of which the case for a 3-4 x 1O 6 M MBH has been clinched, mostly through 
ground-based IR observations of the fast orbital motions of a few stars (Ghez et al., 2003; Schodel et al., 
2003). The second case is NGC4258, which passes a central Keplerian gaseous disc with H2O MASER 
strong sources allowing high resolution VLBI observations down to 0. 16 pc of the centre (Miyoshi et al., 
1995; Herrnstein et al., 1999; Moran et al., 1999). 

In any case, it is nowadays largely accepted that the central dark object required to explain kinematics 
data in local active and non-active galaxies is an MBH. The large number of galaxies surveyed has 
allowed to study the demographics of the MBHs and, in particular, look for correlations with properties 
of the host galaxy. The most remarkable ones are the fact that the MBH has a mass which is roughly 
about 0.1% of the stellar mass of the spheroidal component of the galaxy and that the mass of the BH, 
JZ„ correlates even more tightly with the velocity of this component. These facts certainly strike a 
close link between the formation of the galaxy and the massive object harboured at its centre. 
Of particular importance are the following channels for interactions between stars and the DCO (as- 
sumed to be a MBH). 

The centre-most part of a galaxy, its nucleus consist of a cluster of a few 10 7 to a few 10 8 stars sur- 
rounding the DCO, with a size of a few pc. The nucleus is naturally expected to play a major role in the 
interaction between the DCO and the host galaxy. In the nucleus, stellar densities in excess of 10 6 pc -3 
and relative velocities of order a few 100 to a few 1000 kms~' are reached. In these exceptional con- 
ditions and unlike anywhere else in the bulk of the galaxy, collisional effects come into play. These 
include 2-body relaxation, i.e. mutual gravitational deflections, and genuine contact collisions between 
stars. 

Stars can produce gas to be accreted on to the MBH, through normal stellar evolution, collisions or 
disruptions of stars by the strong central tidal field. When the massive central object acts on the gravity 
of a body, a star, the tidal forces set in and the difference of gravitational forces can prodigally vary 
between the diametrically separated points of the star. The result is that the star will be altered in its 
shape, from its initial approximately spherical architecture to an ellipsoidal one, splitted into two lobes, 
since the volume remains the same. In the end, as tidal forces increase, the star will be tidally disrupted. 
See section 4.13. In Fig. (1.6) we give an intuitive image of this scenario, where distortions due to 
gravitational-lens have not been taken into consideration. In Fig. (1.7), on the left we show a Chandra 
x-ray image of J1242-1 1 with a scale of 40 arcsec on a side. This Fig. pinpoints one of the most extreme 
variability events ever detected in a galaxy. One plausible explanation for the extreme brightness of the 
ROSAT source could be accretion of stars on to a super-massive black hole. On the right we have its 
optical companion piece, obtained with the 1.5 m Danish telescope at ESO/La Silla. The right circle 
indicates the position of the Chandra source in the centre of the brighter galaxy. 

These processes may contribute significantly to the mass of the MBH (Murphy et al., 1991; Freitag and 
Benz, 2002). Tidal disruptions trigger phases of bright accretion that may reveal the presence of a MBH 
in an otherwise quiescent, possibly very distant, galaxy (Hills, 1975; Gezari et al., 2003). 
On the other hand, stars can be swallowed whole if they are kicked directly through the horizon or are 
captured by emission of gravitational waves (GWs). The later process is one of the main targets of the 
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Fig. 1.6: Schematic representation of the tidal disruption process. As the star approaches the central BH, 
tidal forces act on it and tear it apart. Illustration credit M. Weiss 




Fig. 1.7: Optical and x-ray images of RX J1242-11. Credits: ESO/MPE/S.Komossa (left) and 
NASA/CXC/MPE/S.Komossa et al. (right) 



future space-borne GW antenna, LISA (Laser Interferometer Space Antenna). 

For a spherical nucleus in dynamical equilibrium, only collisional effects can bring stars on to the 
"loss-cone", i.e. the very elongated orbits which allow close interaction between a star and the DCO 
(Amaro-Seoane and Spurzem, 2001). 

So far, galactic nuclei have been modelled as isolated spherical clusters (i.e. Murphy et al. (1991); 
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Fig. 1.8: Extension of galactic correlations to smaller systems 



Freitag and Benz (2002)). However, non-spherical structures such as triaxial bulges, bars or stellar discs 
are common on scales of 100-1 000 pc, and also the nucleus itself may be non-spherical. E.g. it could 
be rotating, as a result of a merger with another nucleus (Milosavljevic and Merrit, 2001) or due to 
dissipative interactions between the stars and a dense accretion disc (Rauch, 1995). The influence of 
non-sphericity at small and intermediate scales on the structure and evolution of the nucleus has been 
little explored but it could boost the estimates of capture and disruption rates by orders of magnitudes 
Merritt and Poon, M. Y. (2003). 

1.3.2 Intermediate mass black holes 

If one extends the relations we mentioned above to smaller stellar systems, one could expect that glob- 
ular clusters host so-called intermediate-mass black holes, i.e. BHs whose mass is in the range of 
1O 2 -1O 4 M (see, for illustration, Fig. 1.8). 

After having been suggested in the 70's to explain the x-ray sources observed in globular clusters, 
later discovered to be stellar-mass binaries, this possibility has recently be revived by two lines of 
observations. First IMBHs may explain the ultra-luminous x-ray sources (ULXs) that are present in 
regions of strong stellar formation in interacting galaxies and hence suggesting a link with young "super 
stellar clusters" (SSC), although ULXs are typically not found at the centre of luminous SSCs. On 
IMBHs and their possible link to ULXs, see the review by Miller and Colbert (2003). Second, recent 
HST observations of the stellar kinematics at the centre of M15 around the Milky Way and Gl around 
M3 1 have been interpreted as indications of the presence of an IMBH in both clusters (van der Marel 
et al., 2002; Gerssen et al., 2002, 2003; Gebhardt et al., 2002). However, in the case of M15, the mass 
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of the point masses required by the observations is compatible with zero and TV-body models have been 
made of both clusters that lack a central IMBH but are compatible with the observations (Baumgardt 
et al., 2003a,b) . We note that scenarios have been proposed that would quite naturally explain the 
formation of an IMBH at the centre of a stellar cluster, through run-away stellar collisions, provided 
that the relaxation time is short enough and that very massive stars (1O 2 M < < 1O 4 M ) evolve into 
IMBH (Ebisuzaki et al., 2001; Portegies Zwart and McMillan, 2002; Rasio et al., 2003). 
The theoretical study of the structure and evolution of a stellar cluster (galactic nucleus or globular 
cluster) harbouring a central MBH started 30 years ago. However, due to the complex nature of the 
problem which includes many physical processes and span a huge range of time and length scales, our 
understanding of such systems is still incomplete and, probably, subjected to revision. As in many fields 
of astrophysics, analytical computations can only been applied to highly idealised situations and only a 
very limited variety of numerical methods have been developed so far that can tackle this problem. 

1.4 Super-massive stars: Possible SMBHs progenitors? 

Since it became clear relatively early that most super-massive black holes cannot be formed fast enough 
from stellar mass seed black holes in nuclei (Duncan and Shapiro 1982, 1983, but see also Lee 1995 for 
a somewhat differing view), they must have been formed during the galaxy formation process directly, 
which is linked to cosmological boundary conditions. Rees (1984) argued that galactic nuclei in their 
formation process inevitably produce a dense core consisting of a star-gas system or a cluster of compact 
stellar evolution remnants, both ultimately collapsing to a super-massive black hole. 
The concept of central super-massive stars (SMSs henceforth) {Jl > 5 x 1O 4 M ) 7 embedded in dense 
stellar systems was suggested as a possible explanation for high- energy emissions phenomena occurring 
in AGN and quasars (Hara, 1978; Vil'Koviskii, 1978), such as X-ray emissions (Bahcall and Ostriker, 
1975). SMSs and super-massive black holes (SMBHs) are two possibilities to explain the nature of 
these SMOs, and the super-massive stars may be an intermediate step towards the formation of these 
(Rees, 1984): Stoner and Ptak (1984) argued that a 1O 6 M central gas cloud may be a possible source 
of energy emission of some Seyfert galaxies. For a more up-to-date documentation on the formation of 
such an SMS see Quinlan and Shapiro (1990) (and references therein). On the other hand, such a system 
is interesting not just because it generalises the BH accretion problem to a massive central gas cloud, but 
also because such clouds can be the massive BH progenitors or a possible stock-point to the gas being 
piled up at the centres of GN. Nowadays we know that post-Newtonian instabilities make the lifetime 
of the star lay in the range 10 s — 10 6 yr (Fuller et al., 1986). On the other hand, the active phenomena 
related to QSOs and giant radio sources have a lifetime of about ~ 10 8 yr. This limit caused the idea 
that accretion on super-massive BH could be the engine of these phenomena (Hills, 1975; Lynden-Bell, 
1969; Salpeter, 1964; Lynden-Bell and Rees, 1971b). 

7 Very massive stars are usually defined in the related literature to be stars so massive that at some point they collapse on an 
electron-positron pair instability; according to Bond et al. (1984); Woosley and Weaver (1982); Zeldovich and Novikov (1971), 
this is the situation of a star whose mass is larger than ~ IOOMq but less than 5 x 10 4 Mq. Stars even more massive will collapse 
on the general relativistic gravitational instability before they start burning H (Chandrasekhar, 1964; Hoyle and Fowler, 1963); 
here we denominate such objects after them as super-massive stars, SMSs. 
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Finally, the stability of compact dense star clusters was examined (Zel'Dovich and Podurets, 1966; 
Quinlan and Shapiro, 1990). All these papers as a common feature conclude that, provided the central 
object, star cluster, SMS or a mixture of both, becomes smaller than a certain critical radius, it is able to 
undergo catastrophic collapse in a dynamical time scale due to an instability caused by Post-Newtonian 
relativistic corrections of hydrostatic equilibrium. The question, however, whether and how that final 
unstable state can be reached, is much less clear. Angular momentum of the protogalaxy or its dark 
matter halo, self-enrichment during the dissipative collapse providing opacity through lines which pre- 
vents collapse as compared to radiation driven expansion, and star-gas interactions heating the central 
massive gas object could all at least for some time prevent the ultimate collapse. Given the complex 
physical nature of the interstellar matter, star formation, and stellar interactions alone this is a compli- 
cated question and the conditions under which a super-massive object can form in a spherical, isolated 
star-forming and collapsing gas cloud, rotating or not, has to our knowledge never been exhaustively 
studied and answered. Some pioneering approaches were however done by Spitzer and Saslaw (1966); 
Spitzer and Stone (1967); Colgate (1967); Langbein et al. (1990); Quinlan and Shapiro (1987); Sanders 
(1970). The question has gained even more complexity, since we now know that the baryonic matter 
of galaxies collapses in their dominating dark halo, and that most galaxies and their dark haloes expe- 
rience merging with other dark haloes and large and small galaxies during the hierarchical gravitational 
structure formation (Diaferio et al., 1999; Kauffmann et al., 1999a,b). 

The general idea of how these SMSs have formed is that of massive clusters collapse and coalescence 
nearby the galactic centre. According to Begelman and Rees (1978), a dense cluster of 1M© main 
sequence stars forms a cloud of about 10 5 — 1O 6 M . When the post-Newtonian instability occurs we 
can have either a nuclear explosion, a collapse to BH, fragmentation or a collapse to a BH "seed". 
(Rees, 1984) studies the possible runaway evolution for active GN and includes a scheme in which the 
SMS formation and fate are given (see Fig. 1.9). Fuller et al. (1986) find for their SMSs sample that 
the collapse, bounce, expansion and/or explosion are homologous. They attribute this to the fact that 
the adiabatic index is so close to 4/3 (and the polytropic index n = 3). Goldreich and Weber (1980) 
claim that such a configuration is scale invariant and that it is essentially right at the Jeans mass so that 
collapse or expansion will remain self-similar. 

This has led to another type of study of black hole statistics. As we explained in the last section, the 
black hole masses are well correlated with the bulge masses of their mother galaxies. As a matter of 
fact, these correlations support the idea that black hole formation is linked to galaxy formation. 
The detailed physics and parameters describing how these processes work in a self-consistent model of 
black hole formation, however, are much less understood. This means that we, poor mortals, lack any 
idea of what are the signatures of the black hole formation process in the morphology and kinematics 
of the innermost core and cusp regions, and to what extent they survive the merging history. Brave 
attempts to advance modelling in that domain (Rauch, 1999) demonstrate in our view more the problems 
which still prevail originating from the large dynamical range of the problem and the complexity of the 
treatment of relaxation in a stellar system, rather than that they provide much reliable new insight. 
Galaxy merging poses another serious problem due to the possibility that it can lead to two or more 
black holes in one nucleus, and the structure and kinematics is critically dependent on the evolution and 
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Fig. 1.9: Schematic diagram for the SMS formation and fate (source Rees 1984) 
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possible gravitational radiation merger of the resulting black hole binary. The most direct approach to 
study the dynamical evolution of a system containing a large number of stars is the so-called "direct" N- 
body methods in which the trajectories of N particles are explicitly integrated by computing all N(N — 
l)/2 forces at each time step. With present-day hardware and software, one is limited to N considerably 
smaller than 10 6 and rescaling the results to larger number of stars (10 7 — 10 8 for a galactic nucleus) 
presents intractable difficulties as soon as many different physical processes are at play (relaxation, 
collision, interaction with the central MBH). 

This is the only available possibility here using special GRAPE supercomputers (Makino, 1997; Makino 
and Ebisuzaki, 1996). The use of a GRAPE computer for such simulations is particularly well suited 
because, as a built-in feature, the list of the 50 closest neighbours is returned for each particle, which 
permits an easy local density estimation Burkert (2000). 

The principle behind GRAPE systems is to hardwire on a special purpose chip the most time consuming 
part of an jV-body simulation : the calculation of the accelerations between the particles. The remainder 
is calculated on a normal computer which serves as host to the accelerator board(s) containing the 
special purpose chips. The history of the GRAPE project (short for GRAvity PipE) is detailed in a book 
by Makino and Taiji (1998), and more information is available on the GRAPE project website 8 at the 
University of Tokyo. 

Another possibility are general purpose supercomputers, or a suitable hybrid method between direct and 
approximate A^-body codes (Hemsendorf et al., 2001). From these models it is yet unclear how fast in 
the real system dynamical friction, stochastic three-body interactions and external perturbations work 
together to produce eventually a single black hole again. 

On the other hand, due to the ever increasing observational capabilities with ground and space based 
telescopes we get more and more detailed dynamical and photometric data of the structure of stellar 
systems around black holes. Therefore, we find it worthwhile to reconsider with present day numerical 
possibilities and increased knowledge about galaxy formation and evolution a detailed study of the 
evolution of dense star clusters, with gas, forming an SMS and its further evolution. 
Furthermore, independently on the importance of the role of an SMS in the process of formation of 
SMBHs, these objects have recently drawn the general attention, for they could be good candidates 
for detection by proposed space-based gravitational wave detectors like LISA (the Laser Interferometer 
Space Antenna). LISA would be very sensitive to long wavelengths and low frequency radiation. Thus, 
SMSs are among the most probable sources (see, for instance Thorne 1995). 



1.5 Time-scales 

We introduce in this subsection some useful time-scales to which we will refer often throughout this 
thesis; namely the crossing time, the relaxation time and the dynamical time. These three time-scales 
allow us to delimit our physical system. 

8 http : //grape . c . u-tokyo . ac . jp/grape 
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1.5.1 The relaxation time 

Chandrasekhar (1942) defined a time-scale which stems from the 2-body small-angle encounters and 
gives us a typical time for the evolution of a stellar system. 

This relaxation time could be regarded as an analogy of the shock time of the gas dynamics theory, by 
telling us when a particle (a star) has forgotten its initial conditions or, expressed in another way, when 
the local thermodynamical equilibrium has been reached. Then we can roughly say that the most general 
idea is that this is the time over which the star "forgets" its initial orbit due to the series of gravitational 
tugs caused by the passing-by stars. After a relaxation time the system has lost all information about the 
initial orbits of all the stars. This means that the encounters alter the star orbit from that one it would 
have followed if the distribution of matter were smooth. Therewith can we regard the relaxation time as 
the time interval required for the velocity distribution to reach the Maxwell-Boltzmann form. 
In galactic nuclei the relaxation time is (Chandrasekhar, 1942; Larson, 1970), 

9 a 3 

frelax ~ 160fG 2 mpln(7AT)- (L1) 

In this expression G is the gravitational constant, p the mean stellar mass density, a the local ID velocity 
dispersion, N = \%n c r\ the total particle number and y a parameter of order unity whose exact value 
cannot be defined easily and depends on the initial model and the anisotropy. We will elaborate on this 
in the next chapters. 



1.5.2 The crossing time 

As the name suggests, this is the required time for a star to pass through the system, to cross it. Obvi- 
ously, its value is given by the ratio between space and velocity, 

R 

Across 7 (1-2) 

V 

where R is the radius of the physical system and v the velocity of the star crossing it. 
For instance, in a star cluster it would be: 

fcross — (1-3) 
Oh 

where r^ is the radius containing 50 % of the total mass and Oh is a typical velocity taken at r^. One 
denominates it velocity dispersion and is introduced by the statistical concept of root mean square (RMS) 
dispersion; the variance a 2 gives us a measure for the dispersion, or scatter, of the measurements within 
the statistical population, which in our case is the star sample, 

1 N 

JV i=l 
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Where x, are the individual stellar velocities and pi a is the arithmetic mean, 

1 N 

If Virial equilibrium prevails, we have C7h ~ \/GM^Jy\, then we get the dynamical time-scale 




(1.4) 



where p* is the mean stellar density. 

On the contrary to the gas dynamics, the thermodynamical equilibrium time-scale f re i ax in a stellar 
system is large compared with the crossing time f C ross- In a homogeneous, infinite stellar system, we 
expect in the limit t —* °° some kind of stationary state to be established. The decisive feature for such a 
Virial equilibrium is how quick a perturbation of the system will be smoothed down. 
The dynamical time in Virial equilibrium is (cf e.g. Spitzer, 1987): 

log(yA0 

f dyn x — freiax < freiax ■ (1-5) 

If we have perturbations in the system because of the heat conduction, star accretion on the BH, etc. a 
new Virial equilibrium will be established within a t& yn , which is short. This means that we get again 
a Virial-type equilibrium in a short time. This situation can be considered not far from a Virial-type 
equilibrium. We say that the system changes in a quasi-stationary way. 



1.5.3 Collision time 

fcoii is defined as the mean time which has passed when the number of stars within a volume V = 
L • v re i • Af is one (see Fig. 1.10), where v re i is the relative velocity at infinity of two colliding stars. 
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Computed for an average distance of closest approach r m ; n = |r*, this time is 

«*V(f co ii) = 1 = n*£v re i?coii- (1-6) 

And so, 

Wi= — ^ 1 (1-7) 

with 

E = 7r4 n 1 + — — f ; (1.8) 



mm I 1 _ 2 
\ '"min<-> rel 



C7 r g[ = 2(7^ is the stellar velocity dispersion and E a collisional cross section with gravitational focusing 



1.6 Intention of this thesis 

This work comprises several aspects of theoretical stellar dynamics in clusters, whether globular clusters 
or galactic nuclei, both analytically and numerically. The relevance of computer simulations in current 
astrophysics is well-known and clear, for it is the only way we have to establish (at least in tendency) 
what we should expect from an analytic point of view. Also, questions for which we are not able to 
develop an exact mathematical model are tackled thanks to these methods. Astronomy is a science that, 
at the beginning, was purely observational. We have rarely the chance to access direct measurements of 
what we would like to analyse. In our discipline, creating models that enable us to evaluate the physical 
system exposed to study -although in a way that is more or less idealised according to the method to 
which we resort- is crucial. The work presented here is consecrated to two different aspects that involve 
massive black holes in globular clusters and galactic nuclei. 

As we have seen, the formation of most super-massive black holes (SMBHs) cannot be explained from 
accretion of stars on to seed black holes in nuclei, since the process is unlikely to be fast enough. It has 
been argued that galactic nuclei in their formation process produce a dense core consisting of a star-gas 
system or a cluster of compact stellar remnants. In both cases the system may ultimately collapse to a 
SMBH. Dense super-massive star-gas composite objects have been regarded as a transient progenitor of 
a SMBH. Provided that the central object (a star cluster, SMS or a mixture of both) becomes smaller than 
a certain critical radius, it will collapse in a dynamical time scale due to an instability caused by Post- 
Newtonian relativistic corrections of hydrostatic equilibrium. However, it is not clear at all whether and 
how this final unstable state can be reached. Due to the complex physical nature of the interstellar matter, 
star formation and stellar interactions, this is a rather complicated question. The conditions under which 
a supermassive object can form in a spherical, isolated, star-forming and collapsing gas cloud -rotating 
or not- has never been exhaustively studied. In parallel with my work on stellar clusters containing 
an SMBH, together with an intense analytical study and description of such an SMS (chapter 3), I 
worked on its numerical implementation into the anisotropic model. For that aim, I take into account 
the transfer of radiation in a spherically symmetric moving medium allowing for the contributions which 
are of the order of the flow velocity divided by the velocity of light, the thermal energy equation and the 
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turbulent energy equation. The interaction between stars (collisions) and star gas (drag force between 
stars and gas, loss-cone etc) and post-Newtonian corrections are also considered. How these SMSs could 
power the quasar activity by star accretion and energy flows is one of the main questions that could be 
answered thanks to this method that is in a very advanced state of development, even though not yet 
fully implemented into the numerical program. Here we give an exhaustive analytical description of 
such a physical configuration in chapter 2. I have made also a semi-analytical study of the influence of 
accretion of stars on to an SMS harboured at the centre of the galactic nucleus (see chapter 4). 
For the analysis of star clusters with a central BH, I employed an anisotropic model that solves numer- 
ically moment equations of the full Fokker-Planck equation with Boltzmann-Vlasov terms on the left- 
and interaction (collisional) terms on the right-hand side of the equations. The cluster is modelled like a 
self-gravitating, conducting gas sphere. In this method, all quantities of interest are accessible as smooth 
functions of the radius and time. In chapter 2 there is a detailed description of the approach. 
This enables a detailed study of clean-cut aspects of the dynamics without any noise that particle meth- 
ods suffer from. This model allows us to study the most important physical processes that are present in 
the evolution of a spherical cluster, like self-gravity, two-body relaxation, stellar evolution, collisions, 
binary stars etc and, undoubtedly, the interaction with a central BH and the role of a mass spectrum. I 
have performed calculations to follow the joint evolution of a spherical star cluster with a central BH 
making feasible anisotropy in order to check for the reliability of the method. I include here the study 
of the growth of the central BH due to star accretion at its tidal disruption radius thanks to a diffusion 
model to treat loss-cone physics. The core collapse is studied in detail in a self-consistent manner, as 
well as the post-collapse evolution of the surrounding stellar cluster. The results are in good agreement 
with classical literature about this subject (chapter 5). The current new version of the program enables 
the analysis of the effects of an a (discretised) stellar mass spectrum and stellar evolution. Using it, I 
present also more realistic models of dense clusters and give a description of mass segregation in these 
systems with and without a central BH (chapter 6). 
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Chapter 2 



The theoretical model 



2.1 Introduction 

FUNDAMENTALS and principles of the anisotropic gaseous model are presented in this chapter 1 
so as to give the description of the numerical approach employed to model stellar clusters, as we 
shall see in next chapters. We start with a brief description of the mathematical basis which 
can be regarded as a short summary of the description of kinetic theory in chapter eight of Binney and 
Tremaine (1987). The reader interested in an in-depth study is addressed to it. On the other hand, we 
stress on the peculiarities of the approximation we make use of, the so-called local approximation, so as 
to be aware of the restrictions of our theory and emphasise them. Later, we will give the set of equations 
describing the system, the interaction terms to be taken into consideration to represent realistic models, 
like mass exchange by mass loss, heating, loss-cone, exchange of kinetic energy between the stellar 
system and the interstellar gas via drag forces, etc, to be included in the right-hand terms of this set. 
Ensuingly we shall make a thoughtfully depiction of the gaseous component in the cluster and their own 
interaction terms. 

2.2 The Fokker-Planck equation 

The state of a system of N particles with velocities v = (v;) and positions x = (x ( ) (;' = 1, . . .N) is 
represented by a point in 6A^-dimensional T-space of positional and velocity variables. The jV-particle 

1 At the beginnig of each chapter I will give references of published work related to the contents of the chapter. Typically, I will 
firstly mention the refereed article (if any) and, secondly, contributions to conferences (proceedings). As for the chapter under 
discussion, some sections were used for Amaro-Seoane and Spurzem (2001), Amaro-Seoane and Spurzem (2004) and Spurzem 
etal. (2003) 
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distribution f( N \x,v,t) is defined as the probability to find the system in the volume element d\d\ 
around x and v at time t. is normalised as 1 = / f^d\d\ The evolution of the system is a 

trajectory in T-space; the evolution of a continuous subspace of systems at t = to (initial conditions) can 
be seen as a flow in T-space, which due to the absence of any dissipation is incompressible, and thus 
described by a continuity equation in T-space, which is Liouville's equation: 



dt +L 



N 

I 

i=\ 



IdXi 



f 



(AO 



dxj 
dt 



+ 



,(JV)^ 

1 dt 



0. 



(2.1) 



Using the definition dxj/dt = v, and dvi/dx; = 0, since v, and x, are independent coordinates, and if 
the forces are conservative, dVi/dt = —<?<£>, /<?x,-, where 4>, is the potential at the position of particle ;' 
due to the other particles, we can simplify last equation to 



+E 

i=\ 



dxi dxi d\i 



= 0. 



(2.2) 



Here it has been utilised that dVi/dt does not depend on v, itself, since the potential only depends on the 
spatial coordinates. 

Now we introduce the one-particle distribution function /W (x, v,f) as 



/< 1 )(x,v,f) = J f lN Hx il v i ,t)d 3 x 2 ...d 3 x N d 3 v 2 ---d 3 y N , 
the two-particle distribution function 

/( 2 )(xi,x 2 ,vi.v 2 ,0 = J f {N) (xi,Vi,t)d i x 3 ...d\ N d\ 3 ...d\ N 
and the two-particle correlation function g by 

g(xi,X 2 ,Vl.V 2 ,f) =/ (2) (xi,X 2 ,Vl.V2,f)-/ (1) (xi,Vl,f)/ (1) (x 2 ,V2,f) 



(2.3) 



(2.4) 



(2.5) 



g measures the excess probability of finding a particle at xi, vi due to the presence of another particle 
at X2, V2- Since is normalised to unity, one has p(x2,f) = mN J f^ l \x2,\2,t)d\2, where p is a 
mean mass density, and m the individual stellar mass. Assuming that is symmetric with respect to 
exchange of particles (i.e. all particles are indistinguishable), and observing that 4>i is also symmetric 
with respect to exchanges of the particles 2, . . . ,N, one arrives at 

dfW | ^ dfW N -i dt> { dfW 



dt 



N dxi dv\ 
d 



—Gm(N—l) |_(g(x 1 ,X2,v 1 .V2,f))^-(|x 1 -X2r 1 )A 2 rf 3 V2. (2.6) 
Now one substitutes by the more common phase space density / = /N and drops for simplicity 
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all subscripts "1". It follows 



df d£_N-]_d^_ dj_ 
dt dx N dx d\ 



Sf 
St 



(2.7) 



If the average particle distance d — l/n 1 / 3 is bigger than the impact parameter related to a 90° 
deflection, most of the scattering is due to small angle encounters, which change velocity and position 
of the particle only weakly. So, if we assume that all correlations stem from gravitational two-body 
scatterings of particles, not from higher order correlations, we have that the Fokker-Planck equation is 



Sf 
St 



3 

E 



dx. 



-(/(x,v)D(Ax,))+^-(/(x,v)D(Av i )) 



(2.8) 



Here the convenient notation of diffusion coefficients has been introduced, which contain the integration 
over the velocity and position changes, as e.g.: 



D(Avi) := j Av,- x I / (x,v,Ax,Av)d/ 3 Axd/ 3 Av, 



(2.9) 



where ^(x, v, Ax, Av) d 3 Axd 3 A\dt is defined as the probability for a star with position x and velocity v 
to be scattered into a new phase space volume element ii 3 Ax <i 3 Av located around x + Ax, v + Av during 
the time interval dt. 



2.3 The local approximation 

There are two alternative ways of further simplification. One is the orbit average, which uses that 
any distribution function, being a steady state solution of the collisionless Boltzmann equation, can be 
expressed as a function of the constants of motion of an individual particle (Jeans' theorem). For the 
sake of simplicity, it is assumed that all orbits in the system are regular, as it is the case for example in a 
spherically symmetric potential; thus the distribution function / now only depends on maximally three 
independent integrals of motion (strong Jeans' theorem). Let us transform the Fokker-Planck equation 
to a new set of variables, which comprise the constants of motion instead of the velocities v, . Since in 
a spherically symmetric system the distribution only depends on energy and the modulus of the angular 
momentum vector J, the number of independent coordinates in this example can be reduced from six to 
two, and all terms in the transformed equation (8) containing derivatives to other variables than E and J 
vanish (in particular those containing derivatives to the spatial coordinates x,). Integrating the remaining 



31 



The theoretical model 



parts of the Fokker-Planck equation over the spatial coordinates is called orbit averaging, because in 
our present example (a spherical system) it would be an integration over accessible coordinate space for 
given E and J (which is a spherical shell between r m i n (E,J) and r max (E,J), the minimum and maximum 
radius for stars with energy E and angular momentum J). Such volume integration is, since / does 
not depend any more on x, carried over to the diffusion coefficients D, which become orbit-averaged 
diffusion coefficients. 

Orbit-averaged Fokker-Planck models treat very well the diffusion of orbits according to the changes of 
their constants of motion, taking into account the potential and the orbital structure of the system in a 
self-consistent way. However, they are not free of any problems or approximations. They require checks 
and tests, for example by comparisons with other methods, like the one described in the following. 
We treat relaxation like the addition of a big non-correlated number of two-body encounters. Close 
encounters are rare and thus we admit that each encounter produces a very small deflection angle. 
Thence, relaxation can be regarded as a diffusion process 2 . 

A typical two-body encounter in a large stellar system takes place in a volume whose linear dimen- 
sions are small compared to other typical radii of the system (total system dimension, or scaling radii of 
changes in density or velocity dispersion). Consequently, it is assumed that an encounter only changes 
the velocity, not the position of a particle. Thenceforth, encounters do not produce any changes Ax, so 
all related terms in the Fokker-Planck equation vanish. However, the local approximation goes even 
further and assumes that the entire cumulative effect of all encounters on a test particle can approxi- 
mately be calculated as if the particle were surrounded by a very big homogeneous system with the 
local distribution function (density, velocity dispersions) everywhere. We are left with a Fokker-Planck 
equation containing only derivatives with respect to the velocity variables, but still depending on the 
spatial coordinates (a local Fokker-Planck equation). 

In practical astrophysical applications, the diffusion coefficients occurring in the Fokker-Planck equation 
are not directly calculated, containing the probability *P for a velocity change Av from an initial velocity 
v. Since D(Av,), and D(Av,-Av/) are of the dimension velocity (change) per time unit, and squared 
velocity (change) per time unit, respectively, one calculates such velocity changes in a more direct way, 
considering a test star moving in a homogeneous sea of field stars. Let the test star have a velocity v and 
consider an encounter with a field star of velocity Vf. The result of the encounter (i.e. velocity changes 
Av,- of the test star) is completely determined by the impact parameter p and the relative velocity at 
infinity v re i = |v — V/|; thus by an integration of the type 



the rate of change of the test star velocity due to encounters with v re i, in field of stars with particle density 
n f, averaged over all relevant impact parameters is computed. The integration is normally carried out 
from po (impact parameter for 90" deflection) until R, which is some maximum linear dimension of the 
system under consideration. Such integration generates in subsequent equations the so-called Coulomb 

2 Anyhow, it has been argued that rare deflections with a large angle may play a important role in the vicinity of a BH (Lin and 
Tremaine, 1980). 




(2.10) 
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logarithm In A; we will argue later that it can be well approximated by ln(0.1 IN), where N is the total 
particle number. The diffusion coefficient finally is 

D(Av,) = J <Avi> p f(v f )d 3 v f , (2.11) 

where /(v/) is the velocity distribution of the field stars. In an equal mass system, f(\ f) should be equal 
to the distribution function of the test stars occurring in the Fokker-Planck equation for self-consistency. 
In case of a multi-mass system, however, /(v/) could be different from the test-star distribution, if the 
diffusion coefficient arising from encounters between two different species of stars is to be calculated. 
The diffusion coefficients are (for an exact procedure see Appendix 8. A of Binney and Tremaine 1987): 

D(Avi)= AnG 2 m f \\\A^-h{\) 

D{AviVj)= AnG 2 m f \nA-j-j-g{M) (2.12) 
with the Rosenbluth potentials Rosenbluth et al. (1957) 

h(y)= {m + m f )J J^±Ld\ f 

g(y)= m f J f{v f )\v-v f \d\ f . (2.13) 

With these results we can finally write down the local Fokker-Planck equation in its standard form for 
the Cartesian coordinate system of the vf. 



\StJ enc } [h dvi V >d vj 2 .^j dvtdvj V dvidvj ) 



(2.14) 



Note that in Rosenbluth et al. (1957) the above equation is given in a co variant notation, which allows 
for a straightforward transformation into other curvilinear coordinate systems. 

Before going ahead the question is raised, why such approximation can be reasonable, regarding the 
long-range gravitational force, and the impossibility to shield gravitational forces as in the case of 
Coulomb forces in a plasma by opposite charges. The key is that logarithmic intervals in impact pa- 
rameter p contribute equally to the mean square velocity change of a test particle, provided p^S> po (see 
e.g. Spitzer 1987, chapter 2.1). Imagine that on one hand side the lower limit of impact parameters (po, 
the 90° deflection angle impact parameter) is small compared to the mean interparticle distance d. Let 
on the other hand side D be a typical radius connected with a change in density or velocity dispersions 
(e.g. the scale height in a disc of a galaxy), and R be the maximum total dimension of the system. Just 
to be specific let us assume D = 100c/, and R = 100D. In that case the volume of the spherical shell with 
radius between D and R is 10 6 times larger than the volume of the shell defined by the radii d and D. 
Nevertheless the contribution of both shells to diffusion coefficients or the relaxation time is approxi- 
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mately equal. This is a heuristic illustration why the local approximation is not so bad; the reason is with 
other words that there are a lot more encounters with particles in the outer, larger shell, but the effect 
is exactly compensated by the larger deflection angle for encounters happening with particles from the 
inner shell. If we are in the core or in the plane of a galactic disc the density would fall off further out, 
so the actual error will be smaller than outlined in the above example. By the same reasoning one can 
see, however, that the local approximation for a particle in a low-density region, which suffers from 
relaxation by a nearby density concentration, is prone to failure. 

These rough handy examples should illustrate that under certain conditions the local approximation is 
not a priori bad. On the other hand, it is obvious from our above arguments, that if we are interested 
in relaxation effects on particles in a low-density environment, whose orbit occasionally passes distant, 
high-density regions, the local approximation could be completely wrong. One might think here for 
example of stars on radially elongated orbits in the halo of globular clusters or of stars, globular clusters, 
or other objects as massive black holes, on spherical orbits in the galactic halo, passing the galactic disc. 
In these situations an orbit-averaged treatment seems much more appropriate. 



2.4 A numerical anisotropic model 

In this section we introduce the fundamentals of the numerical method we use to model our system. We 
give a brief description of the mathematical basis of it and the physical idea behind it. The system is 
treated as a continuum, which is only adequate for a large number of stars and in well populated regions 
of the phase space. We consider here spherical symmetry and single-mass stars. We handle relaxation 
in the Fokker-Planck approximation, i.e. like a diffusive process determined by local conditions. We 
make also use of the hydrodynamical approximation; that is to say, only local moments of the velocity 
dispersion are considered, not the full orbital structure. In particular, the effect of the two-body relax- 
ation can be modelled by a local heat flux equation with an appropriately tailored conductivity. Neither 
binaries nor stellar evolution are included at the presented work. As for the hypothesis concerning the 
BH, see section (5.2). 

For our description we use polar coordinates, r 9, <j). The vector v = (v,-) , i = r, , <j) denotes the velocity 
in a local Cartesian coordinate system at the spatial point r,6,(j). For succinctness, we shall employ the 
notation u = v r , v = vg, w = v,),. The distribution function /, is a function of r, t, u, v 2 + w 2 only due to 
spherical symmetry, and is normalised according to 



Here p(r,t) is the mass density; if m* denotes the stellar mass, we get the particle density n = p/m*. 
The Euler-Lagrange equations of motion corresponding to the Lagrange function 




(2.15) 




(2.16) 
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are the following 



d<t> v 2 +w 2 

u = — ^ — I 

or r 

uv w 2 

v= + — -r (2.17) 

r rtan 8 

uw vw 



r ria.nO 

And so we get a complete local Fokker-Planck equation, 

^L +v ^L + v ^L +V - ^L + v - (^-) (2i8) 

dt Vf dr Vf dv v V& dvo V dvq> ySt J pp 

In our model we do not solve the equation directly; we use a so-called momenta process. The momenta 
of the velocity distribution function / are defined as follows 



<i,j,k>:= j vy g \fyf(r,Vi,ve,V4,t)dvidvedvf, 
We define now the following moments of the velocity distribution function, 

< 0,0,0 >:= p = j fdudvdw 

< 1,0,0 >:= u = J ufdudvdw 

< 2,0,0 >:= p v + pu 2 = / u 2 fdudvdw 



(2.19) 



< 0,2,0 >:= p e = j v 2 fdudvdw 

< 0,0,2 >:= ptj, = j w 2 fdudvdw 

< 3,0,0 >:= F r + 3up r + u 3 = J u 3 fdudvdw 

< 1,2,0 >:= Fg+upg = J uv 2 fdudvdw 

< 1,0,2 >:= F^,+up l j ) = J uw 2 fdudvdw, 



(2.20) 



where p is the density of stars, u is the bulk velocity, v r and v t are the radial and tangential flux velocities, 
p v and p t are the radial and tangential pressures, F v is the radial and F t the tangential kinetic energy flux 
(Louis and Spurzem, 1991). Note that the definitions of pi and F, are such that they are proportional to 
the random motion of the stars. Due to spherical symmetry, we have pg = p<j, =: p t and Fg=F^=: Ft/2. 
By p r = pa 2 and p t = pa 2 the random velocity dispersions are given, which are closely related to 
observable properties in stellar clusters. 
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F = (F r +F t )/2 is a radial flux of random kinetic energy. In the notion of gas dynamics it is just an 
energy flux. Whereas for the 0— and <j>— components in the set of Eqs. (2.20) are equal in spherical 
symmetry, for the r and t- quantities this is not true. In stellar clusters the relaxation time is larger 
than the dynamical time and so any possible difference between p r and p t may survive many dynamical 
times. We shall denote such differences anisotropy. Let us define the following velocities of energy 
transport: 

Vr =^~+u, (2.21) 
3 p t 

Ft , 
v t =— + u. 
2p t 

In case of weak isotropy (p T =p t ) 2F V = 3F t , and thus v r = v t , i.e. the (radial) transport velocities of radial 
and tangential random kinetic energy are equal. 

The Fokker-Planck equation (2.18) is multiplicated with various powers of the velocity components u, 
v, w. We get so up to second order a set of moment equations: A mass equation, a continuity equation, 
an Euler equation (force) and radial and tangential energy equations. The system of equations is closed 
by a phenomenological heat flux equation for the flux of radial and tangential RMS (root mean square) 
kinetic energy, both in radial direction. The concept is physically similar to that of Lynden-Bell and 
Eggleton (1980). The set of equations is 



dt r 2 dr ' 



du du GM T 1 dp r t - Pi- Pt „ 

— h u -z — I H r 1-2 = 

dt dr r z p dr pr 



dpr Id. 9 du 1 d , 9„. 

37 + -j-( r 2 upi)+2 Pt — + - — {r 2 F r ) (2.22) 
at r L or or r L or 

2F t _ 4 (2p r -p t ) 

r 5 A^frelax 



-di + ^Tr {r UPl)+2 — + 2r^J-r {rFl) 
Ft = 2 (2p r - p t ) 

r 5 Xa t relax 

where Xa is a numerical constant related to the time-scale of collisional anisotropy decay. The value 
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chosen for it has been discussed in comparison with direct simulations performed with the A^-body code 
(Giersz and Spurzem, 1994). The authors find that Xa =0.1 is the physically realistic value inside the 
half-mass radius for all cases of N, provided that close encounters and binary activity do not carry out 
an important role in the system, what is, on the other hand, inherent to systems with a big number of 
particles, as this is. 

With the definition of the mass M r contained in a sphere of radius r 

?=Wp, (2.23) 
or 

the set of Eqs. (2.22) is equivalent to gas-dynamical equations coupled with the equation of Poisson. To 
close it we need an independent relation, for moment equations of order n contain moments of order 
n + 1. For this intent we use the heat conduction closure, a phenomenological approach obtained in 
an analogous way to gas dynamics. It was used for the first time by Lynden-Bell and Eggleton (1980) 
but restricted to isotropy. In this approximation one assumes that heat transport is proportional to the 
temperature gradient, 

dT da 2 

F = -K— = -A^— (2.24) 
or or 

That is the reason why such models are usually also called conducting gas sphere models. 
It has been argued that for the classical approach A oc X 2 /%, one has to choose the Jeans' length 
Xj = <J 2 1 \4nGp) and the standard Chandrasekhar local relaxation time f re i a x <*■ C 3 /P (Lynden-Bell 
and Eggleton, 1980), where X is the mean free path and T the collisional time. In this context we obtain 
a conductivity A oc p/o. We shall consider this as a working hypothesis. For the anisotropic model 
we use a mean velocity dispersion a 2 = {a 2 + 2a; 2 ) /3 for the temperature gradient and assume v r = v t 
(Bettwieser and Spurzem, 1986). Forasmuch as, the equations we need to close our model are 



X do 2 

\% Gp f r elax dr 



= 
= v t . 



(2.25) 



2.5 Interaction terms and implementation of the gaseous compo- 
nent 

The set of Eqs. (2.22) describe the stellar cluster without taking into consideration the different inter- 
action processes which may play a role in these systems. Heating, cooling processes, mass exchange 
between stars and gas by stellar mass loss or disruptive stellar collisions, exchange of kinetic internal 
energy between the stellar system and the interstellar gas via drag forces, loss-cone processes etc are 
not included in the equations. Here we will give a description for some interesting phenomena that crop 
up in the systems we want to study. Chapter 5 is partly devoted to a profound description of a loss-cone 



37 



The theoretical model 



diffusion model. 

In parallel with the description of the interaction terms for the stellar component of a cluster, here we 
will describe the inclusion of an SMSin the anisotropic numerical model. For that aim, one should 
take into account the transfer of radiation in a spherically symmetric moving medium allowing for 
the contributions which are of the order of the flow velocity divided by the velocity of light, the thermal 
energy equation and the turbulent energy equation. The equations will be given also "in their logarithmic 
form" (i.e. in our equations we work with lnx instead of x), ready to include in the gaseous model. 
The interaction terms for the gaseous component 3 are also described and post-Newtonian corrections 
are also considered. How these SMSscould power the quasar activity by star accretion and energy flows 
is one of the main questions that could be answered thanks to this method. Therefore we have found of 
paramount importance to give the set of equations describing this situation. 



2.5.1 The star component 

We now introduce the interaction terms to be added to right hand of the star component equations. This 
and the next section benefit from the paper by Langbein et al. (1990). 



Equation of continuity 

Langbein et al. (1990) derive the interaction terms to be added to the basic equations of the gaseous 
model. According to them, the star continuity equation is no longer 



dp* 1 d 2 
dt r 2 «?r 



-(r 2 p*u*)=0, (2.26) 
ct r- ox 

but 



■ + 3T(r P*u+) = -=- + "ST : ( 2 - 27 ) 



at r^ V y coll \5t 

where the right-hand term reflects the time variation of the star's density due to stars interactions (i.e. due 
to the calculation of the mean rate of gas production by stars collisions) and loss-cone (stars plunging 
onto the central object). 

If /(vrei) is the stellar distribution of relative velocities, then the mean rate of gas production by stellar 
collisions is 

' 8 -£) =-f ^^/(v rel )^ 3 v rel (2.28) 

01 / coll • / kel[>°~coll 'coll 



3 One should here not be confused by the terminology. When we talk about the gaseous model, we mean the numerical 
anisotropic model we use to model stellar systems; it is called like that because it is based on moments of the Boltzmann equation 
with relaxation. The stellar cluster is modelled like a self-gravitating, conducting gas sphere, and was first presented in Louis 
and Spurzem (1991) and Giersz and Spurzem (1994). On the other hand, when we refer to the gaseous component, we mean the 
physical gas in the cluster, which may conform an SMS. 
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In the calculation of equation (14) /(v re i) is a Schwarzschild-Boltzmann distribution, 

(v re l,r - «*) 2 



/(Vrel) 



1 



•exp 



(2.29) 



As regards f c , it is the relative fraction of mass liberated per stellar collision into the gaseous medium. 
Under certain assumptions given in the initial work of Spitzer and Saslaw (1966), we can calculate it as 
an average over all impact parameters resulting in r m ; n < 2r* and as a function of the relative velocity at 
infinity of the two colliding stars, v re i. Langbein et al. (1990) approximate their result by 



/c(v re l) = 

with <7 co ii = 100. So, we have that 

The first interaction term is 
St 



( 1 + <?coll \/ Ocoll / Vrel) 1 V re l > O co ll 
V re l < Ocoll, 



/c(Vrel) = 



0.01 Ocoll = V rc l 







Ocoll > Vrel, 



fc 



coll tcoll 

which, for simplification, we re-call like this 

St 



1 - erf 



( Ocoll 

VV6cJ r 



1-erf 



( Ocoll 

\V6ot 



-P*X C0 11. 



(2.30) 



(2.31) 



coll 



Since the evolution of the system that we are studying can be regarded as stationary, we introduce for 
each equation the "logarithmic variables" in order to study the evolution at long-term. In the other 
hand, if the system happens to have quick changes, we should use the "non-logarithmic" version of 
the equations. For this reason we will write at the end of each subsection the equation in terms of the 
logarithmic variables. 

In the case of the equation of continuity, we develop it and divide it by p* because we are looking for 
the logarithm of the stars density, d\np+/dt = (l/p*)dp*/dt. The result is: 



o\np± du± 
dt dr 



d\np* lu* 
or r 



p* v 5t 



+ W8_p ± 

coll P* V St 



(2.32) 



Momentum balance 

The second equation in Eq. (2.22) has the following star interaction terms: 



du* du+ GMr 1 dp r p T - p t 
• + U* — 1 ^ 1 r r Z- 



dt 



P* dr 



( 8u± 
\~St~ 



(2.33) 
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The interaction term is due to the decelerating force at which stars that move inside the gas are subject 
to. As we shall see, an estimate for the force is given by Eq. (4.6). Explicitely, it is 



8uA 1 



= -Xdrag — - Kg) (2.34) 



where we have introduced the following definition: 



^drag = -C — -p*P g <7tot, (2.35) 



with a t Q t = a r 2 + a 2 + (m* — M g 



To the end of the calculation of the logarithmic variable version of the equation, we multiply Eq. (2.33) 
by p±r/p r : 

p±r f du+ \ GM r d\np r p t , 



, ■>■+«* + p*+ 3T ^+2(l-^) = -X drag -(K*-w g ) (2.36) 

p r \ dt J rp r dlnr p r p x 



Radial energy equation 

As regards the last but one equation, the interaction terms are: 
where 



= -2X drag cr r 2 , f - -X coll p,d r 2 e. (2.38) 

5t J drag V Of J coll 

In order to determine £ we introduce the ratio k of kinetic energy of the remaining mass after the 
encounter over its initial value (before the encounter); k is a measure of the inelasticity of the collision: 
for k — 1 we have the minimal inelasticity, just the kinetic energy of the liberated mass fraction is 
dissipated, whereas if k < 1 a surplus amount of stellar kinetic energy is dissipated during the collision 
(tidal interactions and excitation of stellar oscillations). If we calculate the energy loss in the stellar 
system per unit volume as a function of k we obtain 

e = /- 1 [l-fe(l-/ c )]. (2.39) 

We divide by p r so that we get the logarithmic variable version of the equation. We make also the 
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following substitution: 



F r = 3p r v r 

F t = 2p t v t (2.40) 



The resulting equation is 



■ + (M* + 3v r ) — h3 — + — +- t/* + 3v r -2v t — + 



dt dr \ dr dr J r \ ' p, 



5 ?relax p r \ 8t J , p r \ dt 



coll 



Tangential energy equation 



To conclude the set of equations of the star component with the interaction terms, we have the following 
equation: 



^ +(%) . (2-42) 



& 7,w V 5f 



coll 



where 

--2X drag (J t 2 , ^ = -X collP *G t 2 e. (2.43) 

We follow the same path like in the last case and so we get the following logarithmic variable equation: 



dlnp t . „ ,d\np t d . „ . 4. „ , 

+ w* + 2v t ) + 3- («* + 2v t ) + - w* + 2v t - (2.44) 

<7f or dr r 

42f t -l _ l f8 Pt \ + ]_fd Pt 



5 frelax Pt \ 5f / drag Pt V 5f /coll 

2.5.2 The gaseous component 

In this section we give the set of equations corresponding to the gaseous component as for their right 
hand interaction terms. 

Equation of continuity 

For the SMS the equation of continuity looks as follows: 
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where, for the mass conservation, we have that, obviously, 

5 Pg\ ( $P* 



&t I coll V ^ f / coll 



(2.46) 



We follow the same procedure as for the star continuity equation to get the equation in terms of the 
logarithmic variables: 

ainp g du g ainp g 2u g _ 1 / 5p g \ 

— + ^7 + u *— + — -p- g {-W) con (2 ' 47) 

The interaction term is in this case 

s(£)-s(-£)~£*- 

Momentum balance 

We modify equation number (2.9) of Langbein et al. (1990) in the following way: 

d(p gUg ) dp g du g , 

dt ~ g dt +Pg dt ' ( y) 

we substitute this equality in their equation, divide by p g (u g is the variable in our code) and make use 
of the equation of continuity for the gas component. Thus, we get the following expression: 

+Mg ^ + ^ + - = ( °£) (2.50) 

dt b dr r z Pg or c \ Ot 



coll 



To get the interaction term we use the mass and momentum conservation: 

'8p g \ i ( 8p^ 



St ' coll \ 8t / coU 



= 

n 

8(PgU g )\ , /" 5(p*W*) 



We know that 



thus, 



» ! J coll ' V & 7 coll °" (2 ' 51) 

: 0, (2.52) 



5m* 
"57 



coll 



, = w*p*X co n = p g I ) +M g X co iip*. (2.53) 

5f /coll V « / coll 



Therefore, the resulting interaction term is 



^) =— ^coiiK-Mg) (2.54) 

°' / coll Pg 
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In the case of the stellar system 

F= l -(F, + F t ) = 5 -p^ (2.55) 



By analogy, we now introduce F Ta ^ in this way 



^rad _ _ 5 



H = -PgV g , (2.56) 



where v g is per gas particle. 

vg = ^- (2.57) 

As means to write the equation in its "logarithmic variable version", we multiply the equation by p g r/ p g , 
as we did for the corresponding momentum balance star equation and replace H by 5/2p g v g , 

p„r ( du„ du„\ GM V d\np z 5 fC ex t 

p g \dt h dr J rp g ^ h dlnr 2 c rfe 8 

r 

— p*^coii(«*-Mg) (2.58) 
Pg 

Radiation transfer 

We extend here and improve the work done by Langbein et al. (1990) by resorting to a more detailed 
description of the radiation transfer (Castor, 1972). 

Consider a radiation field; we place a surface element da with a surface normal n, see Fig. (2.1); the 
radiation energy which passes through da per unit time at angle 9 to n within a small range of solid 
angle d(0 given by the directional angles 9 and <j) is 

dE =I v (d,<l))ndvdad(0 (2.59) 

where jU = cos 9. 

The radiation intensity I v (d,<j>) is defined as the amount of energy that passes through a surface normal 
to the direction (9,<j>) per unit solid angle (1 steradian) and unit frequency range (1 Hz) in one second. 
The intensity of the total radiation is given by integrating over all frequencies, 



/o 



/'OO 

/ I v dv (2.60) 

Jo 

The three radiation moments (the moments of order zero, one and two) are defined by: 

J= \ J v dv = / dv- / I v d[l 

Jo Jo 2 J -i 



1 



+ i 



H= / H v dv= dv- I v ixdn (2.61) 

Jo Jo 2J_i 

K= I K v dv = / dv- / I v ll 2 djl, 

Jo Jo 2 J -i 
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Fig. 2.1: Radiation intensity definition 



The moment of order zero is related to the density of energy of the field of radiation E a ^, the moment 
of order one to the flux of radiation F ra d and the moment of order one to the radiation pressure /? rac i, 





A% 


£rad — 


—J 




c 


^rad = 


AkH 




A ^K 


Pad = 


c 



The transfer of radiation in a spherically symmetric moving medium is considered taking into account 
the contributions which are of the order of the flow velocity divided by the velocity of light; we include 
also the variation from the centre up to the atmosphere of the Eddington factor / E( jd = K/J, where K 
and J are the radiation moments; / E( jd is obtained from a numerical solution of the equation of radiation 
transfer in spherical geometry Yorke (1980). We get the radiation transfer equations by re-writing the 
frequency-integrated moment equations from Castor (1972): 



IdJ | dH | 2H 7(3/ E dd-l) i; J{l+fEdd)dlnp g = 
c dt dr r cr c dt 

K abs (B-J) 



(2.63) 
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IdH d(Jf Edd ) 7(3/ E dd-l) 2« g 2 glngg 

— ^ — I ^ 1 ti r — ti — 

c at dr r cr c at 

-KenPgH (2.64) 
In the equations K" a b s and JC e xt are the absorption and extinction coefficients per unit mass 

K"abs = 8 J — ~i ""ext = Pgl^abs + Kscatt)) (2.65) 
a 

A(r) is the cooling function, B the Planck function and fc scatt the scattering coefficient per unit mass. 
We have made use of dM r /dr = 4n 2 p, /edd = K/J, and the Kirchhoff's law, B v = j v /K v (j v is the 
emission coefficient), so that the right-hand terms in Castor (1972) are the corresponding given here. 
We now look for the logarithmic variable version of both equations; for this aim, we divide Eq. (2.63) 
by / and multiply Eq. (2.64) by 2c/ (5p s v s ), 

1 dlnJ 5 d 5 3/ Edd - 1 

;( v gPg) + 7ZPs v z ~ u g-{ l +fiw 



dt ' 2Jdr K ***' Jr** * 

(B-J) (2.66) 



1 dlnp g _ K" abS/ 



c dt J 



dlnv g | d\np g | 2c 1 d(Jf E dd) | 2c 3/ E dd-l / 2t< g 
dt dt 5 p g v g dr 5 rp s v g r 

~ dln Pg 



dt 



-cfQxtPg (2.67) 



Where we have substituted H = 5p g v g /2. 



Thermal energy conservation 

It is enlightening to construct from Eqs. (2.22) an equation for the energy per volume unit e = (p r + 
2p t )/2 which, in the case of an isotropic gas (p r = p t ) is e = 3p/2. For this aim we take, for instance, 
equation (2.37) and in the term 2p v du+/dr we include now a source for radiation pressure, 2(p T + 
Pvnd)du+/dr and we divide everything by e so that we get the logarithmic variables. The resulting 
equation is 

dine . „ dine 2 f . 1 f 8e\ I / 8e\ ,„ „ N 

-^— + (u g + 3v g )^— + -(u g + v g ) = -( — \ +- — 2.68) 
dt s g; dr r x g g; e\8t) Am e\8t J coll 

The interaction terms for this equation are 

(w) d = W<*r 2 + Oi 2 + («*-«g) 2 ) (2-69) 
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Fig. 2.2: Representation of the logarithmic radial mesh used in the code. With v, p we represent both, the 
radial and tangential components of the velocity and pressure. 



5e\ 



^W>*((a r + a t ) 2 e + (k* - M g ) 2 - E, a 2 oll ) 



(2.70) 



Mass conservation 

The mass conservation is guaranteed by 



3 dMr 



471 d. 



- = P* + Pg 



(2.71) 



2.6 A mathematical view of our approach 

In this subsection, we explain briefly how the gaseous model is solved numerically. We concentrate on 
aspects of the method not exposed in previous work. This description is therefore complementary to 
Sec. 2.2 of Giersz and Spurzem (1994). The algorithm used is a partially implicit Newton-Raphson- 
Henyey iterative scheme (Henyey et al. 1959, see also Kippenhahn and Weigert 1994, Sec. 1 1.2). 
Putting aside the bounding conditions, the set of equations to be solved are Eqs. 2.22 to 2.25. In practice, 
however, the equations are rewritten using the logarithm of all positive quantities as dependant functions. 
As explained in Giersz and Spurzem (1994), this greatly improves energy conservation. Formally, one 
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may write this system as follows 



dt 




for i = 1 ... 4 



fori = 5...JV e 



cq 



(2.72) 



where the *W are the local quantities defining the state of the cluster, i.e. 



(2.73) 



= {logp,M,log/? r ,log/? t ,logM r ,v r -t/,v t -t/} 



with N cq = 7 in the present application. 

To be solved numerically, this set of coupled partial differential equations have to be discretised accord- 
ing to time and radius. Let us first consider time stepping. Let At be the time step. Assume we know the 
solution x(t — At ) at time t — At and want to compute x(t). For the sake of numerical stability, a partially 
implicit scheme is used. We adopt the shorthand notations *W = x^(t) and = x^(t — At). Time 
derivation is replaced by finite differences, 



In the terms /W, we replace the x^ by which are values intermediate between and x^\ x^ = 
fr^ + (1 - QyU\ with C = 0.55 for stability purpose (Giersz and Spurzem, 1994). 
Spatial discretisation is done by defining all quantities (at a given time) on a radial mesh, {r\ , ri , . . . r^ r } 
with r\ = and r^ r = r max . A staggered mesh is implemented. While values of r, u, v t , v r and M r 
are defined at the boundaries of the mesh cells, p, p t and p r are defined at the centre of each cell, 
see Fig. 2.2. When the value of a "boundary" quantity is needed at the centre of a cell, or vice-versa, 
one resort to simple averaging, i.e. bt = 0.5(/3j ; _i + b^), c# = 0.5 (c* + c^+i), if b and c are is border- 
and centre-defined quantities, and b, c their centre- and border-interpolations, respectively. For all runs 
presented here, N T = 300 and the points r2 to r max are logarithmically equidistant with r max = 10 4 pc 
and r2 — 1.7 x 10~ 6 pc. Let us adopt the notation xj/' for the value of x^ at position r# (or and 
Ark = — r^-i . Then, radial derivatives in the terms /W are approximated by finite differences, 



dxV 



dt 



At 



(2.74) 



(2.75) 



dr Ark 



47 



The theoretical model 



if the derivative has to be evaluated at a point where Xk is defined (centre or border of a cell), or 



dj) _ dj) 

x k x k- 1 

An 



2An 



(2.76) 



otherwise. As an exception we use upstream differencing in du/dr for the second equation in set 2.22, 
i.e. the difference quotient is displaced by half a mesh point upstream to improve stability. 
By making the substitutions for dx^/dt and dx^/dr in the set of differential equations 2.72, one 
obtains, at each mesh point r k , a set of N eq non-linear algebraic equations linking the new values to be 
determined, x Jc _ 1 and x k , to the "old" ones, y k l and y^, which are known, 



^ t (0 (a-i,s l |y t _ 1 .y i )=o 

i = l...N eq , k=l...N t . 



(2.77) 



Note that the structure of the equations is the same at all mesh points, except k = 1 and k = N T . In 
particular, terms k — 1 do not appear in . Also, one has to keep in mind that only the x k _ x and x^ 
are unknown; the y kx and y k play the role of fixed parameters in these equation (as do the Ar k ). If one 

defines a (N eq x iV r ) -dimension vector JT* whose component N eq (k — 1) + ;' is x k \ one can write the 
system of A^ eq x N r equations as — 0, i.e. 
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(2.78) 



The system is solved iteratively using Newton-Raphson scheme. If i£"*[ m ] is the approximation to the 
solution of Eq. 2.78 after iteration m, with [ m j = ( JT* [ m j ) ^ 0, the solution is refined through the 
relation 



^ \m+\\ — 



(2.79) 



48 



2.6 A mathematical view of our approach 



where (d&* /d SE*) 1 is the inverse of the matrix of derivatives. The latter, of dimension (N eq N T ) x 
(N e qN r ), has the following structure 



I ■ □+ 

□ _ ■ D-i 

□ _ 



V 



□ 



□_ 



□_ 



(2.80) 



In this diagram, each square is a N eq x Neq sub-matrix. For 2 < k < N T — 1, lines N eq k — 6 to N eq k of 
d,^* jd X* are composed of a group of 3 such N eq x A^q matrices, that span columns 



N eq k — 13 to N eq k + N eq , while the rest is composed of zeros, 
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(2.81) 



We can see this more explicitely, 
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More schematically, to have an overview of the structure, 




The Heyney method is a way to take advantage of the special structure of matrix d.jp* /d3£* to solve 
system (2.79) efficiently, with a number of operation scaling like 0(N r ) rather than ff{N^) as would 
be the case if one uses a general-purpose matrix inversion scheme 4 . Setting = — &*{ m \ and W* = 

4 Memory usage is also reduced, scaling like (?(N r ) rather than ff(N^). 



50 



2.6 A mathematical view of our approach 



&*[m+\\ ~ 3£*[m], Eq. (2.79) is equivalent to 



(2.82) 



with W* the unknown vector. We further decompose vectors W* and 3§* into jV eq -dimensional sub- 
vectors, each one representing the values at a given mesh points, 

r 2 



n 

Then, the system (2.83) can be written as a set of coupled jV eq -dimensional vector equations, 



(2.83) 



(2.84) 



The algorithm operates in two passes. First, going from k = 1 to N r , one defines recursively a sequence 
of iVeq-vectors % and (N eq x Af eq )-matrices 9Jl k through 



m l = (m l )- 1 D +l 

m k = (m k -n- k im k - l )- 1 n +k 2<k<N r . 



(2.85) 



(2.86) 



TIn, is not defined. In the second pass, the values of the unknown % are computed, climbing back from 
k = N T to 1, with 

n = %-mn+r \<k<N r -\. 

Note that, with this algorithm, only (N &q x Af eq ) matrices have to be inverted. We use Gauss elimination 
for this purpose because this venerable technique proves to be robust enough to properly deal with the 
kind of badly conditioned matrices that often appear in this application. 

The initial model for the Newton-Raphson algorithm is given by the structure of the cluster at the 
previous time, 2£* m (t ) = JT* (t — At) One iterates until the following convergence criteria are met. Let 
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us set 8xfi = x$ 



[m+l] 



,K ! J " r *=l..JV r 

with ei = 10~ 6 . For velocities (m, v r — h, v t — u), one checks 

i ^ / 5# ^ 



. Then, the condition for logarithmic quantities is 

i=l..JVeaiVr V * / 



max - V — — * < e 2 , (2.88) 

i=l...iVe, A^r jjf^ \ X M + £lWfc i 



v 4' +£iW kj 

with £2 = 10~ 3 and Wk = r^AnGpk) 1 ^ 2 ■ Generally, two iterations are sufficient to reach convergence. 
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Chapter 3 



Super-massive stars 



3.1 On the nature and peculiarities of a super-massive gaseous ob- 
ject 

ARISING from 1962, we find in the literature one of the first proposals on this topic 1 : Some 10 8 
supernova? were thought to happen in 10 6 years; since they were also thought to be the nor- 
mal galactic source of cosmic rays, we would have found the origin for this strong emission. 
Nonetheless, this was a too strong requirement for Burbidge and Burbidge (1962), who found the argu- 
ment not convincing. He proposed the idea that one supernova in a tightly packed set of neighbouring 
stars could trigger explosions in these. If Burbidge and Burbidge (1962) felt satisfied with such an idea, 
F. Hoyle and W. A. Fowler could not say the same, for they found it not convincing. To fulfil Burbidge's 
argument the stars should be so close packed together that they would be physically in contact. It was 
impossible for F. Hoyle and W. A. Fowler to withstand the seduction of the simple idea that at the centres 
of galaxies star-like objects exist with masses of up to 1O 8 M (Hoyle and Fowler, 1963). As regards the 
stability and/or existence of such an object, they "turn a blind eye". Actually, we have to wait one year 
to find a deep study on this subject by Chandrasekhar. 

F. Hoyle and W. A. Fowler developed their argument with the argument that wholly convective stars 
could "do the job": Consider an ideal gas with radiation pressure: 

P = P gas +P r =^pr + V, (3.1) 

'Part of this chapter was used for Just and Amaro-Seoane (2001), Amaro-Seoane et al. (2002) and Just and Amaro-Seoane 
(2003) 
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P, jj,, p, and T being total pressure, mean molecular weight, density and temperature; 3i = k/m u — 
8.315 x 10 7 erg/ (Kg) is the universal gas constant, k = 1.38 x 10~ 16 eig/K the Boltzmann constant, 
m u = 1 amu = 1.66 x 10~ 24 g the atomic gas unit and a = 7.565 x 10~ 15 erg/(cm 3 Zf 4 ) the radiation- 
density constant. We use here the gas constant with a dimension (energy per K and per unit mass); in 
thermodynamics we usually find that it is energy per K and per mole. Because of that, the molecular 
weight ji is dimensionless (instead of [jit] = mass/mole). Assume that j3 = Pgas/P is constant throughout 
the whole star; thus, we can see in 1 — j3 = aT 4 / (3P) that j3 = constant implies the relation T 4 <~ P. 
Therefore, 

which is a polytropic relation with a polytropic index n = 3 for constant j3 . The polytropic constant 

is a free parameter since we can choose j3 in the interval [0,1]. 
Eddington's quartic equation is 

The left-hand side of the equality becomes smaller in the interval [0, 1] when M increases. In other 
words, the radiation pressure becomes more important the larger the stellar mass is. This means that 
SMSs are dominated by the radiation pressure. As a consequence, the adiabatic gradient V ac j is reduced, 
(particularly V ac j — > 1/4 when ^3 > 0, see e.g. Kippenhahn and Weigert 1994) and the star becomes 
convective with V = V a( j. This adiabatic structure requires constant specific entropy s. Thanks to the 
first law of thermodynamics, we have that 

4aT 3 

s=-r-. (3.5) 
3p 

Constant specific entropy means p ^T 3 and from the pressure equation we can see that P ~ p ; therefore 
P ~ p 4 / 3 , which again means that SMSs are polytropes with n = 3. 

The SMSs are star polytropes with a free K, which implies that their mass M can be chosen arbitrarily. 
For each mass, 1 — J3 /(jU 4 j3 4 ) can be obtained and then, Eq. (3.3) gives us the corresponding value of 
K. But if the mass is given, there still exist an infinite number of models for different radius (in spite 
of the fact that K is already determined by M) since K is independent of R. With a specific entropy 
being constant, we have that the adiabatic polytropic constant y is roughly the local adiabatic index 
T=[d{\nP)/d{\np)] s . 

The SMS evolves radiating energy and entropy and, in the case of contemplating rotation, losses also 
angular momentum via mass shedding. Since the pressure is dominated by radiation, the luminosity of 
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an SMS is close to the Eddington limit, 



% = ^Edd 



= 3.4x 10' 



,4 




(3.6) 



0Th 



where Jz? is the luminosity of the SMS , m p the proton mass and a T h the Thomson cross section. 

3.1.1 Nuclear energy source 

The most important nuclear energy source for SMSs is hydrogen burning. We can distinguish between 
two regimes: For temperatures less than T = 5 x 10 & K, hydrogen burns on the standard j3 -limited 
CNO cycle (Fowler, 1965). For temperatures greater than the above-mentioned one, Wallace and 
Woosley (1981) proved that leak out of the j3 -limited CNO cycle can be very important, because of 
the 15 0(a, y) 19 Ne reaction. When we are out of the cycle, we have a flow consisting of radiative proton 
captures and positron decays which builds up toward the iron peak: this is the so-called rp-process. 
Wallace and Woosley (1981) proved that when this process operates, it has 200 or 300 times the energy 
generation rate of the j3-limited CNO cycle. Nevertheless, for masses in excess of a few 1O 5 M , nuclear 
burning is irrelevant prior to the gravitational instability when we compare the nuclear energy genera- 
tion rate with the photon luminosity (see e.g. Fowler 1966; Zeldovich and Novikov 1971), and electron- 
positron annihilation are not important (New and Shapiro, 2001); the creation of electron-positron pairs 
can be neglected forM > 1O 4 M (i.e., T < T ait < 2.5 x 10 9 K) (Bond et al., 1984). 

3.1.2 Instabilities of radiation-dominated stars 

The bulk of the SMSs is expected to be convective, thus the entropy is roughly constant. For stars in 
which the specific entropy is nearly constant, the adiabatic polytropic constant is approximately the local 
adiabatic index (Chandrasekhar, 1939). The evolution of these objects is a rather complicated subject, 
for these stars have always a pressure averaged adiabatic index of by 4/3. In order to keep hydrostatic 
equilibrium, the central temperature must increase as the star mass rises up. The ratio of gas pressure P g 
to the total pressure P = P v + P gas (where P r is the radiation pressure) for an index n = 3 polytrope is 



where ^# is the mass of the SMS. For stars with a huge mass like these objects, the gas pressure is a 
minor perturbation on the total pressure of the star, P g can only be a minute fraction of the total pressure. 
The condition for hydrostatic equilibrium can be expressed as a variational principle (Chandrasekhar, 
1964; Hoyle and Fowler, 1963); the total energy of the star, which is a function of entropy and 
central density is an extremum. We can thus get the equilibrium mass for a given entropy and central 
density by extremising this total energy. If we employ a low enough central density, the equilibrium 
energy will be almost zero, because of the y ^ 4/3. The SMS will quasi-statically contract, radiating 
away entropy and increasing the central density while the equilibrium energy decreases. Gravity will 




(3.7) 
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have a minimum due to general relativistic corrections, so that beyond a certain p cr i t , energy must be 
supplied in order to get equilibrium. The SMS radiates quasi-statically its entropy away and shrinks to 
the point where the central density reaches the p cr j t and instability sets in. At this point, the equilibrium 
energy is stationary, and we have instability, for a second derivative of the total energy of the system is 
zero. Thus, the SMS becomes dynamically unstable and collapse sets in. 

The nature of the SMSs is such that small effects that are normally negligible in the study of stellar 
evolution must here be included to determine the stability properties. 

• general relativity 

• electron-positron pair formation 

• rotation 

• dissociation of heavy nuclei 

In the framework of general relativity, dynamical instability happens before a gaseous mass shrinks 
to the limiting radius compatible with hydrostatic equilibrium. A star whose structure is determined 
almost completely by Newtonian gravitation turns out to be unstable because of general relativity. The 
instability in radiation dominated stars stems from the fact that 7 is nearly 4/3; they "are trembling on 
the verge of instability" (Fowler, 1964). This was first suggested by Chandrasekhar (1964) and Feynman 
(1963) and first applied to stars by Hoyle and Fowler (1963). 

From Chandrasekhar (1964), (post-Newtonian) instability sets in when the radius of the star is less 
than a critical radius & cnt . He shows that if the ratio of specific heats 7 = C p /C v exceeds 4/3 only by a 
small amount, then dynamical instability will occur if the mass contract to the radius 



where Jff is a constant depending on the density distribution in the configuration and the term between 
brackets is the Schwarzschild radius of the SMS. Chandrasekhar evaluates the J(f value for the homoge- 
neous sphere of constant energy density and the polytropes of indices n = 1, 2, 3. The ratio of specific 
heats 7 can be shown to be in the limit of high entropy (Chandrasekhar, 1939), 



to first order in B. This means that the radiation-dominated stars have a 7 close to the value 4/3 and 
then they are on the limit of instability. With the last formula and n = 3 we can approximate the & cnt by 




(3.8) 



4 B , 

y«3 + f +^ 2 ), 



(3.9) 




(3.10) 
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The Schwarzschild radius of the SMS is: 



— = 10 [Wm q )^ (3 - n) 

Then, with the j3 formula, we can find out that 

/ \ 3 / 2 

^= 0A ^[mr Q ) (3 - 12) 

On the other hand, the critical central density for a radiation- dominated, high entropy, index n = 3 
polytrope is 

p cnt = 2xl0-(i-) 3 (^) 7/2 g /cm 3 , (3.13) 

(Fuller et al., 1986). For a higher central density, the nonlinear effects of gravity have a destabilising 
influence to radial perturbations. The critical central temperature is 



r crit = 2.5 x 10 13 ( I K, (3.14) 



l 

M . 

and the corresponding equilibrium energy is 

£ crit = -3.6 x 10 54 erg, (3.15) 

which is independent of . 



3.2 Fencing in the existence zone of an SMS 

It is interesting to examine the time-scales of SMSs in order to know whether or not they are on hy- 
drostatic equilibrium, and so to have a general overview of their possible evolution. The time scale 
associated with hydrodynamic processes in radiation-dominated stars is 

^^v^ w2xl0 "(^) 7/4s - (3 ' 16) 

In order to see a hydrostatic evolution phase, this time scale must be shorter than the modified Kelvin- 
Helmholtz time, also called "cooling time" or "contracting time" for it is the time-scale which gives us 
an approximated idea of how long can live the the star from its gravitational energy. In the life of a star 
there are phases in which this is the main or even the only stellar energy source, and so we say that "the 
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star evolves on a Tkh". This time is given by 




3 x 10 



16 



M 



) 



(3.17) 



where Lgdd is the Eddington luminosity and E m [ n is the value of the total energy of the star at the 
relativistic or electron-positron pair instability point (Shapiro and Teukolsky, 1983). We find that 



For instance, for a 10 M© SMS, Td yn ~ 6.33 yrs and Tkh ~ 9.5 yrs. For times which are comparable to 
the time-scale we have to take into account the cooling of the SMS, since it can determine its evolution. 
For short times compared to the time-scale (t -C T K h) we can forget about the cooling when we want 
to study the state and/or the evolution of the SMS. We have to evaluate which process is more effective 
(which process has a shorter time-scale), either the cooling, the accretion of mass onto the central object 
or the heating of the SMS due to un/confined stars. In such a case, the evolution will not be determined 
by the cooling, but by the faster processes. 

In the case that Tkh > %yd, the cooling is slower than the setting of the hydrodynamical equilibrium and 
so we get a quasi- hydrostatic SMS. At this point there are two completely different aspects which should 
be studied in detail. The first of them is the state of equilibrium. We have a hydrostatic equilibrium, 
we consider that we have no cooling and P rac i » P g . The Schwarzschild criterion for dynamical stability 
says that a star with V ra( j > V ac i (the late being = 4/3, n = 3 polytrope) will be unstable, supposing a 
homogeneous chemical composition of the star, so that = and the small perturbations will increase 
until we have a fully convective star. If the equilibrium is stable, we have a contracting SMS, but we are 
close to instability; on the other hand, if it is not stable after the post-Newtonian instability, we have a 
free-fall collapse, because the radiation pressure grows slower than the gravity force when it contracts. 
The second aspect with which we have to deal is the quasistationary evolution, which is dominated now 
by Tkh an d by the slow cooling. Ditto, we have a slow cooling and contracting SMS provided that we 
get such a state of equilibrium after having analysed the foregoing questions. To complete the study we 
have to know whether the stationary state is stable and how can we describe it. 
We can illustrate this situation thanks to the gravitational and total energy (E gr , E tot ), which are 




(3.18) 



Etot — t 




(3.19) 



(3.20) 



In these equations, ^ pc is the radius of the SMS in pc and Jt% its mass in units of 1O 8 M . We have 
also used the fact that 



3.3 Possible origin of a BH: sequence of SMSs 



As we already mentioned, the contraction rate of the SMSis determined by the Eddington luminosity 
which, expressed in (see eq. 3.6), is 

jSf 8 = 1.3 - 10 46 ^g— (3.22) 

s 

Thus, we have a contraction time (Kelvin-Helmholtz time) 

Tkh = ^ = 3.7M 8 1/2 J R pt 1 yr (3.23) 

L Ed 

To determine whether a SMS can exist or not (from a point of view of the hydrostatic equilibrium), we 
just have to check for its radius in a plane. We can define the radius 3%% at the place where Tkh is 

equal to the free-fall time. In fig. (3.1) the time-scale Tkh m which an SMSevolves is shown, from right 
to left. This gives us an up-limit to the limits in the plane for the SMS. As for the low-limit, we have to 
take into account the condition imposed by the post-Newtonian unstability M w or go even "lower" and 
look for lower masses when nuclear burning at the centre sets in, 2%^. The area where such an object 
can exist (is stable) starts at Ma .When hydrogen burning at the centre sets in at M B b, or when Post- 
Newtonian unstabilities triumph over thermal pressure stabilisation, at & pn , the SMSbecomes unstable. 
For masses < 5 • 1O 4 M , the electron-positron pair instability occurs. The friction time-scale is also 
plotted (Tf r ic = 10 5 , 10 3 , 10 1 yr from right to left) for a typical star with v* = 1500 km/s (assuming that 
1 % of the orbit crosses the inner part of the SMS). For radii smaller than the coupling radius ^ couv \ such 
stars are effectively decelerated. Some typical values of Tm c f° r s °l ar type stars are also plotted. 



3.3 Possible origin of a BH: sequence of SMSs 

It may be enlightening to distinguish two competing processes to have a first look of the evolution of 
this physical problem: The contraction (or collapse with re-bounce) of the mixed core with stars that are 
trapped by friction within the gas and slow down additionally to build a new highly condensed stellar 
core. This core may re-heat the gas or decouple to build a new core. A core-halo interaction is also 
possible. The loss-cone stars of the surrounding stellar system heat the core and feed it by means of new 
trapped stars. This core collapse can be conked for a while until the loss-cone is "empty" (not rigorously 
speaking, but for the practice we can assume it is empty) or the core becomes too massive. 
The friction force of the stars moving through the SMS has a dynamical and a geometrical component. 
For SMSs more compact than the coupling radius & couv i, green line in Fig. (3.1), the typical friction time- 
scale is shorter than the contraction time of the SMS. Thus, stars are effectively decelerated by dynamical 
friction. As a result, loss-cone stars will be segregated to the centre building a highly concentrated sub- 
system, which may essentially change the evolution of the SMSby enhanced interactions. 
The result may be a sequence of SMS -stars cores contained one inside the other until the relatively low 
mass innermost SMS initiates hydrogen burning or collapses to a BH "seed". In Fig. (3.2) we give an 
intuitive scheme for this scenario. 



61 



Super-massive stars 




* (pc) 

Fig. 3.1: Characteristic radii and time-scales for an SMS. The Schwarzschild radius £?schw ' s f ar below the 
"allowed" region of an SMS. The evolutionary tracks are horizontal lines from right to left. The contraction 
time Tkh is also displayed in order to show that without re-heating the life time of the SMS is at 2- 10 4 yr 
before nuclear burning or the post-Newtonian collapse sets in (see text) 



3.4 Stabilisation theory 
3.4.1 The role of rotation 

Rotation is likely to play a decisive role in the quasistationary evolution of an SMS, as well as in its final 
collapse. Since SMSs are on the verge of stability, rotation can prolongate their equilibrium evolution. 
The gravitational instability sets in when Tsms yields a value less than a critical value, which is roughly 

rsMS<r ^iz| + 1 , 2 2W (3 , 4) 
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where Mschw = 2G^/c 2 is the Schwarzschild radius of the star and rj = T/\W\ the ratio of the ro- 
tational energy to the gravitational potential energy. For 77 = 0, the last equation becomes r cl -j t = 
4/3 + 1.12^schwA^' (Misner et al. 1973 and references therein). The adiabatic index grows thus as 
the radius of the star becomes smaller. Rotation has a stabilising effect and can hold the collapse if 

„ 1 4-3(r SMS -1.12^ Schw /^) 
n >rj crit ~ 25 _ 3(rsMS _ 1 A2 ^/M) ^ 

Baumgarte and Shapiro (1999b, a) study the effects of rotation on a SMS and follow the evolution of 
these objects up to the dynamical instability and collapse. They find out that ratios like St/Msam, 
St I r/|W| and Jc/{G^£ 2 ) (T is the rotational kinetic energy, W the gravitational potential energy 
and J the angular momentum) for a maximally and rigidly rotating n = 3 polytrope at the onset of radial 
instability are universal numbers, "key ratios", as they denominate them, which are independent of the 
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mass, spin, radius or even the history of the star. They are thus universal constants: 

(wX„r om - G&L-"* 

with the polar radius Sl v w 7M e /3 (8% e stands for the equatorial radius). This deformation could be 
responsible for the reduction of the luminosity by about 36% below the usual Eddington luminosity 
from a non-rotating SMS (Baumgarte and Shapiro, 1999b,a): 

Lrot = 0.639L E dd. (3.27) 



Since luminosity determines the timescale (i.e. the mean life) of these objects, rotation seems to be a 
factor whose importance should be seriously taken into account. 

They perform a fully relativistic, numerical calculation and compare it with an analytical treatment, 
which is in good agreement. They claim that cooling by photon radiation drives the evolution, while 
the rotating SMS loses mass, angular momentum and entropy. Gas pressure may stabilise SMSs in 
the absence of angular momentum (Zeldovich and Novikov, 1971; Shapiro and Teukolsky, 1983), but 
Baumgarte and Shapiro in their paper maintain that even a small degree of rotation could dominate gas 
pressure. They analyse also the final collapse of rotating SMSs, which might not necessary be a BH. 
They reckon that the collapse may be inhibited, in the form of a disc and/or bar, or possible fragmenting 
into several blobs. 

The study of the importance of a dark matter background on the stability of an SMS has also been done, 
with the conclusion that although it has stabilising effects, it can be neglected compared with the effects 
of rotation (Bisnovatyi-Kogan, 1998). 

Nevertheless, all the work being performed as regards the stability of SMSs in the case that such objects 
would rotate, and other topics, like the equilibrium and stability of SMSs in binary systems assume that 
the system is an isolated one; namely, the accretion and heating are assumed to be negligible. 



3.4.2 Stabilisation by fluctuations 

The density, temperature and velocity fluctuations excited by the stellar component in the gas (Just 
et al., 1986; Deiss et al., 1990) do not only decelerate the stars by dynamical friction, but there is also 
a stabilising effect on the gas component, which can be described by an effective j3. The diffusion of 
photons is much faster than the dynamical time scale in the SMS. Since the temperature is determined by 
the radiation field, the perturabations are isothermal as long as the matter is optically thick. For scales 
smaller than the mean free path A = 3.3 10~ 6 M% R^.p c /p pc, the density fluctuations are adiabatic 
leading to an enhanced mean thermal pressure. The basic corrections to the unperturbed parameter j3o 
are 

£ = j3o(l+/^) + j3 turb >j3 (3.28) 
Po 
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where j3 tur b gives the ratio of turbulent to radiation pressure and (p£) is the unperturbed density. The 
scaling constant K and the fluctuations' density (pf) must be determined by some fluctuation theory. 
Using a quasilinear approximation of a mode analysis the perturbations excited by the stars moving 
through the gas will be estimated to determine the magnitude of the corrections due to the turbulent 
velocity field and the adiabatic density fluctuations on the stability parameter j3o- An increase of the 
effective j3 will shift the line for the Post-Newtonian instability to higher masses. 
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Chapter 4 



A semi-analytical approach to dense 
gas- star systems 

4.1 Introduction 

DYNAMICS in galactic nuclei of dense stellar systems is the subject to be presented in this chapter 
1 for some of the physical configurations we introduced in former sections. It has been often 
suggested that, in addition to the thermonuclear sources inside the constituent stars, the active 
phenomena of high-energy emissions occurring at galactic nuclei and quasars could be explained resort- 
ing to the idea that a massive dark object is being harboured in the centre of the galaxies (Magorrian 
et al., 1998). 

We will start this discussion facing the black hole star accretion problem in a dense stellar medium in a 
given time scale. Radii of stars are extremely small (r w 10 10 cm) in comparison with the mean typical 
distances between stars in astrophysical stellar systems (a typical interstellar distance in our galactic 
neighbourhood is w 10 18 cm). Even though stellar collisions occur seldom, a single star will modify 
its orbit at random because of the small pushes it receives from stars which pass close by. We can 
roughly define a certain time required for the stars "to forget their initial orbits" due to two-body, small- 
angle gravitational encounters, the gravitational tugs of these passing-by stars; this would be one of our 
timescales of interest, the relaxation time. 

Then will we follow through the initial problem introducing the concept of super-massive star. As a 
first description of it, a super-massive star > 1O 4 M , where ^ is its mass) is a dense gas cloud 
located in the centre of dense stellar systems which can be regarded as a general case for the black hole 

'Part of the work exposed in this chapter was used for Amaro-Seoane and Spurzem (2001a), Amaro-Seoane and Spurzem 
(2001b) and Amaro-Seoane and Spurzem (2000) 
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accretion problem. 

If the super-massive star scheme embraces the black hole accretion problem, we could think that the 
most general case is to be analysed first, since it supposedly gives an overview of the general problem. 
Nevertheless, if we dare to envisage the possibility of doing it upside-down, if we generalise the example 
and try to extract from it the global concept, it reports us the satisfactory surprise that everything turns 
out to be clearer after having evaluated the particular case. The most general case becomes pointless if 
it yields as a result a more difficult understanding of the subject submitted at study. 

4.2 General concepts 

In this section we discuss the consequences from the astrophysical and dynamical point of view of the 
presence of a black hole (or a massive compact central object, from now onwards just BH) in a dense 
stellar object in about a relaxation time. We consider the steady-state distribution and consumption of 
stars orbiting a massive object at the centre of a spherical, stellar system. 

The distribution of stars is determined by the relaxation processes associated with gravitational stellar 
encounters and by the consumption of low angular momentum stars which pass within a small distance 
of the central mass. Stars whose orbits carry them within the tidal radius r t of the BH will be tidally 
disrupted. The peribarathron 2 (distance of closest approach to the BH) is determined by the specific 
orbital angular momentum L and by the BH mass. There are situations in which the stars that have a 
radially elongated orbit and a low angular momentum pass close by the system centre and interact with 
the massive central object. In such a situation it is interesting to evaluate the density of those stars whose 
angular momenta are limited by a superior L m ; n ; that is, the stars which belong to a defined region in the 
velocity-space. 

Stars at a position r whose velocities are limited by a superior limit vi c (r) and, consequently, with an 
angular momentum L < L m ; n = rv\ c , have orbits that will cross the tidal radius of the central BH in 
their motion. They will be disrupted due to the tidal forces and then they are lost for the stellar system. 
Such stars are said to belong to the loss-cone, since they are lost for the stellar system. The loss-cone is 
depleted in a crossing time f cross = rj a, where a is the ID velocity dispersion. 

The diffusion of stars into this loss-cone has been studied by (Frank and Rees, 1976) and by Lightman 
and Shapiro (1977). We will talk about a "critical radius" within which stars on orbits with r < r cr j t 
diffuse. Inside they are swallowed by the BH, after being scattered to low angular momentum loss-cone 
orbits. 

If there is a central point mass ./#., such that » m*, then its potential well will affect the stellar 
velocity field out to a distance 

r h = G^./<7 2 . (4.1) 

This expression gives us the influence radius of the central object. G is the gravitational constant. 

The star gets disrupted whenever the work exerted over the star by the tidal force exceeds its own binding 

2 This word fits quite well the idea of closest approach to this "sinking" hole for it has the meaning of no return. The barathron 
{P&paQpov) was in the ancient Greece a cliff down to an unreacheable or unseen place where criminals were thrown 
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energy. If we compute the work exerted over the star by the BH we can get an expression for the tidal 
radius (see Addendum A of this chapter, section 4.13), 



n = 



1 V3 



-(5-n) 

3 fn* 



r*, 



(4.2) 



where n is the polytropic index of the star, the mass of the star and its radius. 



4.3 Loss-cone phenomena 

Frank and Rees (1976) studied how a stationary stellar density profile around a massive star accreting 
BH looks like. They found that the density profile follows a power-law within the region where the 
gravity of the massive star dominates the self-gravity of the stars, 



p oc r 1/4 . (4.3) 

This was followed by intensive numerical studies by other authors (Shapiro and Marchant, 1978; Marchant 
and Shapiro, 1979, 1980; Shapiro and Teukolsky, 1985) which are all in agreement with the first work 
of Frank and Rees. 

The replenishment of the loss-cone happens thanks to the small-angle gravitational encounters in a 
timescale which is for the most real stellar systems slower than the dynamical processes. 
We have loss-cone effects also in the neighbourhood of a massive central gas-formed object, an SMS, 
see Fig. (4.1). Stars with such orbits enter the gas-formed central object and lose kinetic energy if 
their density is high enough. Nevertheless, such stars will not disappear from the stellar system just 
by one crossing of the central object, as it happens for the BH problem. In this scenario the stars lose 
their energy in each crossing and their orbits come closer and closer to the central object until they are 
"trapped" in it, their orbits do not extend further than the massive central object radius (confined stars). 
This process was described qualitatively by Hara (1978), da Costa (1981) and Hagio (1986). Da Costa 
suggested the name "dissipation-cone". However, we will keep the loss-cone term for it is the most 
commonly found in the related literature. We have to take into account that the meaning is that of a 
defined region of the velocity space at the position r, even though there is no quick loss. 
In Addendum B of this chapter (section 4. 14) we derive the following expression for the diffusion angle 
(the mean deviation of a star orbit in a dynamical time fd yn ): 

0d - J£±. (4.4) 

V 'relax 

Now we look for a condition at a place r > n, for a star to touch or to cross the influence radius of the 
central object within a crossing time. For this aim we look now for the amount of stars which reaches 
the central influence radius with an unperturbed orbit. Unperturbed here means that the star orbit results 
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from the influence of the gravitational potential and from that of the rest of the stars and of the central 
object, and it is not affected by the local, two-body, small- or big-angle gravitational encounters. We 
envisage then the average part of the gravitational potential, whereas the random component due to the 
individual behaviour of the stars will be neglected. 

To define the loss-cone angle we say that a star belongs to this cone when its distance to the peribarathron 
(which depends on the orbit we have, i.e., on the energy E and angular momentum) is less or equal to 
the tidal radius, r p (E,L) < r t , 9 < 9\ c . 

There is a maximum for which the peribarathron radius is less than or equal to the tidal radius. We 
define this as 9\ c (where the subscript "lc" stands for loss-cone). In Addendum B we derive an expression 
for it in terms of the influence radius of the central object. 

As for the Vi c (r), we get its value using angular momentum and energy conservation arguments in 
Addendum C (section 4.15). 

v lc (r) = ■ 72[<Hr t )-<K'-)]+v r (r)2. (4.5) 

4.4 The critical radius 

It is interesting to evaluate a certain radius which Frank and Rees ( 1 976) introduced by defining the ratio 
4 : = (he/ 0d- When E, = 1, then 9\ Q = 9®, and this corresponds to a "critical radius", r cr j t , if there is only 
one radius with this condition. Inside the critical radius (i.e. E, > 1, 9[ c > 9o) stars are removed on a f^yn- 
For larger radii (i.e. E, < 1, 9\ c < 0d) we cannot talk about a "loss-cone" because this 0d corresponds to 
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the variation of 9 within a fdy n , and this is the required time for the star within the loss-cone to plunge 
onto the BH if the orbit is unperturbed. The angle variation happens sooner than the required time for 
stars to sink into the loss-cone. If 9d> 9i c , loss-cone stars can get in and out of the loss-cone faster than 
they could reach the central object. 

4.5 The loss-cone contribution to the heating rate of the gas 

Dissipation of the stellar kinetic energy of a star plunging onto the SMS and suffering from the drag 
force leads to a heating of the SMS. Another possible consequence of the local star-gas interaction is 
the formation of massive stars within the cloud due to the accretion of ambient gas (da Costa, 1979). 
This could increase the supernova rate and be an important source of energy. Here we will assume an 
equilibrium for the SMS for the timescales of interest. 

The star distribution will be affected at large radii (r » 8%, 8% being the SMS radius) by removal of stars 
in the central regions of a stellar system (Peebles, 1972). A drift of stars occurs in the centre of the stellar 
system in order to recover the equilibrium. In the special case of having a BH at the centre of the system, 
processes like tidal disruption lead to the destruction of the star. In this arena we have an outward energy 
flux created by the sinking of stars via relaxation processes (local, two-body, small-angle gravitational 
encounters). 

When we consider the general case of an SMS, the basic picture is the same, but some aspects vary; the 
removal of stars and the inward transport are due to a different process. The inward flux of energy is 
produced not only for the local, two-body, small-angle gravitational encounters, but also for dissipative 
processes. The effective sinking for stars is now related to gas-drag energy dissipation and can involve 
or not their actual physical destruction. A star moving through the cloud will quickly dissipate energy 
because of gas drag and then it will sink into the centre of the SMS. The main difference between a BH 
and an SMS is basically that the former produces a low-angular momentum star depletion on a crossing 
time, whereas the dense gas cloud or SMS does the same but in a dissipation time. 
Since we want to analyse the effects on the dense stellar system arising from the presence of a central gas 
cloud we have to distinguish between those stars whose orbits are limited to the region where the SMS is 
located (confined stars) and those stars in orbits which surpass the radius of the cloud (unconfined stars). 
When we talk about confined stars, the first steps in the evolution are determined by the energy dissipa- 
tion given by the drag force that the dense gas cloud exerts on the individual stars. Stars lose velocity in 
their motion inside the cloud, they are slowed down by the gas and therefore they cede heat to the cloud. 
The slowing down of the stars makes them become a more compact subsystem which will sink down to 
the centre of the cloud. The system becomes self-gravitating and we have a cusp in the stellar distribu- 
tion. However, star-star interactions can play a decisive role in this point, since they yield a depletion 
in the number of confined stars, or we can have direct collisions between them and thus disruption or 
coalescence. This could avoid that a singularity in the core collapse crops up. We cannot exclude the 
exchange of mass between individual stars and the gas as another possible way to prevent the singularity, 
since this can yield the star disruption via stellar wind, or the creation of heavier stars which become a 
supernova (da Costa, 1979). 
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The consequent evolution of the confined system will be in part determined by the rate at which sur- 
rounding stars outside the gas cloud refill this confined-stars gap. The importance of the core collapse 
will also be a decisive point for the evolution. 

Unconfined stars move on orbits extending larger than the SMS radius. However, they can suffer its 
influence out to a radius within which the presence of the SMS is effective. The idea is exactly the same 
as for the BH, for the SMS is a generalisation of the former case 



4.6 Kinetic energy dissipation 



The drag force that the individual stars suffer when they cross the SMS is given by the next equation 
estimated by Bisnovatyi-Kogan and Syunyaev (1972): 

F D = C D Sp SMS vl (4.6) 

Cd is a numerical parameter of order unity, Psms is the mean density of the gas cloud, v* is the velocity 
of the stars and S is the cross section of the stars, S = %r^. In case of a supersonic motion of the star 
the force Fo can be interpreted as caused by the ram pressure (pressure difference) originating at a bow 
shock in front of the moving star. Due to the physical shock conditions one can show in such a case 
C D w 4 (Courant and Friederichs, 1998). 

Suppose that the star crosses the SMS from one extreme to the opposite, i.e. along its diameter; thus, if 
& is the radius of the SMS, the stellar energy dissipated during each passage through it is 



A£ D = Frj • 2@. (4.7) 

The orbits of the stars within will be elliptic shaped with one focus at the SMS centre. The semi-major 
orbit axis a will shrink because of the drag force, driving the orbit directly into the SMS. The average 
energy dissipation rate is 

dE = A£ D = 2C D psMsnrlG^ 
dt " T Ky/4a*/Gje ' 

where T is the period. 

In the BH accretion problem the individual stars are doomed to sink in the hole if they belong to the 
loss-cone, whereas for the SMS it depends on the number of passages through the gas cloud. It means 
that the non-confined star can cross the gas cloud several times until it is trapped. To make this picture 
clearer we resort to the BH, particularly when we introduced the concept of 9d- In the case of SMS, 

9d „ /2H2L^ZE. (4.9) 

V frelax 

We can get the number of necessary passages for the star to be trapped if we write the potential difference 



74 



4.7 Loss-cone stars velocity field distribution function 



between r and Si (for the BH problem it is n = 1), 

A</> = G^fm*(i-^). (4.10) 

Then, the number of passages is: 

« trap = |A4>|/A£ D . (4.11) 

The dissipation of stellar kinetic energy defines a region in phase-space which we will call the loss-cone 
too, by analogy with the BH accretion problem. We have to recall however, that the star will be blown 
up by tidal forces and converted into gas by just one passage inside the tidal radius r t for the case of 
having a BH whereas the number of passages is n trap > 1 if we consider the gas cloud situation. 



4.7 Loss-cone stars velocity field distribution function 

It stands to reason that at distances much larger than the SMSradius (r » SI) the star field velocity 
distribution has a Maxwellian shape: 

f(r,v)= g(r)-h(v), 

2 2 

f(r,v)= AT(r).exp(-^)ex P (-^) (4.12) 
Taking into account that d 3 v = 2nv t dv r dv t , the density is 

P( r )=J f(r,v)4nv 2 dv = n^N(r)o I 0?; (4.13) 

then 



Pir) 

And so get we the normalised field velocity distribution, 



N(r) = — QJL _ 2 (4.14) 



^■•"■^(i^^^-4?^ < -ik } - <4 ' 15) 

In order to get the density of stars within the loss-cone we have to compute the following integral: 

Pk{r)= / /(v r ,v e ,v ( | ) )dv r (ivec/v0. (4.16) 

J — "r|max ''—Vic 

Taking into account that Avgdv^ = 2nv t dv t and that f(r, v) = f(r, — v), 

/-v r | max /•vi c (r,v r ) 

p k (r) =4np{r) I dv r / /(r,v r ,v t )v t dv t . (4.17) 
Jo Jo 
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In this expression, the maximal radial velocity is v 



r|max — ^escape 

of both, the super-massive star and the stellar system potential, 

0W = 0sms(O + 0*W- 
The integral happens to be analytical, and it yields the following result: 

pic(r,Vr)=p(r)-{a-£-p-\l/}, 

where 



= y/2<j>(r) and the potential is the sum 



(4.18) 



(4.19) 



a = 
P = 
W = 

A0 = 



ed{^{r)/a}) 

j Sf- 2Ad>\ 

erf(i^0(r)/<T r 2 ) 

/ {r 2 -M 2 )a 2 
\{r 2 -@ 2 )o 2 + @ 2 o 2 

0(^)-0«- 



1/2 



(4.20) 



Since we are working with a Gaussian function whose width is a, the contributions of velocities v r > 2a r 
to the total mass are small and therefore negligible. In the practice it means that we can approximate 
the integral by v r = 2a r . In such a situation, the loss-cone adopts an easy geometrical form: an open 
cylinder in the — v r direction with a radius vi c . The resulting integral yields 



Plc(r) 
P(r) 



1 — exp < — 



2a 2 + A(p(r) 



•erf 



'0W 



(4.21) 



We use for the stellar system a polytrope with n = 5, with a compact core and an extended outer envelope 
(named Plummer model after Plummer (1911), since it was first used by him to fit the observed light 
clusters distributions), which is commonly used for several numerical model computations. 
In this model the potential is spherical, 







(4.22) 



where Jl± is the system's total mass and r c its core radius. 
Thus, 



(r 2 + r 2 ) 5 / 2 ' 



(4.23) 
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Then, if we resort to the equation of Poisson, we can find out that this system happens to have the 
following physical properties: 

3M / r 2 \- 5 / 2 
ft»«= 1^3 1 + 3 



c \ c 

r 3 />f 



Thus, 



»■ yjr 2 + r 2 c ' 

Therewith, the resulting expression is 



0W = — — + T - (+ 25 ) 



= (l-expA)-erf(B), (4.26) 

p(0 



where 



A= - 




GJ( GJt 



+ r 2 ^Jr 2 + rl 

The velocity vectors of all stars which belong to a given (fixed) phase space density f = fo = const 
shape an ellipsoid whose two tangential and one radial major axes have lenghts equal to the velocity 
dispersion: Gg = o§, and G r (Frank and Rees, 1976). If one uses fo — /((7 r , Gg, Oq), the surface of the 
ellipsoid A = %GgG^G t is a measure for the available velocity space. 
A cone of angle 9] c , 

9i c :=arcsin (— ), (4.28) 

cuts out a segment of the foregoing velocity ellipsoid's A surface of corresponding fraction surface A\ c , 

A Xc ^%dl, (Sic «%). (4.29) 

We can then define the ratio Q! := A\ c /A rts 0^/4, which is a measure for the loss-cone size. 

With the assumption that vi c does not depend on v r any more and taking into account that with a 

Schwarzschild-Boltzmann distribution we can reduct the loss-cone momenta to elementary Gaussian 
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error functions. The quantity 




(4.30) 



yields in first order 




4 



(4.31) 



This is the connection between the preceding simple picture of the loss-cone and the definition of Q. in 
the velocity space: For a Schwarzschild-Boltzmann distribution we can find out that, at first order, Q. is 
of the same size as £1' . 

4.8 Isotropy and anisotropy in the stellar system 

We introduce now an isotropy ratio in order to study the different possible situations for the stellar 
distribution: The tangential velocity dispersion is O 2 = c? + O^; in case of isotropy, o r 2 = O? + Og, then 



Now we define the ratio R := 2o r 2 / O t 2 . According to this definition, R = 1 for the isotropic case. The 
corresponding values of R for radial and tangential anisotropy can be obtained bearing in mind that 



o 2 = o 2 + o 2 = o 2 (/?/2 + 1), o t = v/y/R/2+1 and o r = o t • y/R/2. 

We have the loss-cone star density as a function of the super-massive star radius, S%. This does not 
provide much information, since in principle this radius could have any size; we do not have a criterion 
for it yet. Instead, what does make sense is to express this loss-cone star density in terms of the super- 
massive star stability, which is something has been studied in detail (Fuller et al., 1986). 
From Chandrasekhar (1964), instability sets in when the radius of the star M is less than a critical radius 
3% znX . He shows that if the ratio of specific heats y = C p /C v exceeds 4/3 only by a small amount, then 
dynamical instability will occur if the mass contract to the radius M cnt 



Thus, we introduce the stability coefficient 8 :=M/M cnt . We just have to substitute M = 8-^ cnt in the 
loss-cone star density formula and vary 8 instead of 

4.9 Connection at the influence radius 

Since we are interested in the diffusion angle, we now derive two expressions for it, within the influence 
radius of the SMS and outside it. For this aim we look at the dynamical and relaxation time at this radius. 

The o(r) varies depending on r^. 






(4.32) 



or<r b : 
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In this case we assume that the gravitational potential is dominated by the SMS. Therefore, the velocity 
dispersion is: 



<rW = y— ; (4-33) 

This is just an approximation which we make here for simplicity, since we stand on radii r < r^. To 
include r<n,we need a better model, which can only be obtained by numerical solution of the equation 
of Poisson, and this is subject of future work. 

or>r b : 

For radii larger than rj, we will use a Plummer model for the velocity dispersion, 



However, we have to match both solutions, within and outside the SMS influence radius, for we have to 
look for a velocity dispersion connection; otherwise we get artificial, non-physical "jumps" in the plots 
for r m r^. This can be performed as follows: We add in the velocity dispersion expression a factor a, 



ff(r) = V^T (i+^)^ ' (4 ' 35) 

and we ask now for both velocities dispersions to be equal at the influence radius: 

a{r b )\ r<rh = a{r h )\ r>rii . (4.36) 
We can conclude that a has the following value: 



" = \& ( ^ ,,/4 <437) 

Note that a is necessary because our velocity dispersion is approximate for r < i\, but not for r <C i\. 
Therefore, the velocity dispersion outside the influence radius is: 



a{r) 



r c rh 



r 2 



1/4 

(4.38) 



As regards the dynamical time, as we saw in the introduction, fdyn = r/ C7 r (r), is 

,■3/2 . 

r <Hi 



IGJ( 

1 .. <h I \ r >r h 



GJ( ' \ r^+rl 



For the relaxation time we need the density of the stellar system; We already commented that, according 
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to Frank and Rees (1976) or to Bahcall and Wolf (1976), the density of the spherical stellar system 
around the central mass should follow a power law of the form p (r) oc r~ 7 / 4 . 

/ \ - 7 / 4 

p(r)=p p (r h ). - (4.40) 



where p p (Vh) is 



Therefore, 



r 2\-5/2 



c \ 7 c 



3-#* I ' . n 



2\ -V2 / r \ -7/4 



P« = 7^3(l + 3) (-) (4-42) 



47^ \ r\ 

Thus, the final expression for the relaxation time in such a system is 

3y^ c 3 /G^\3/2/ r x7/4/ r^5/2 



i 4^1n(yiV) v , , v , h / v / r4 4T> 

— - /^3/2 fl A5/2r % ( 2 2^3/4 ^™ 

These expressions enable us to get the value for 0q; we can get 9\ c thanks to £1, 9[ c w 2\/ii. Then we 
plot them altogether to find out at what radius do they cross, since we are looking for the critical radius, 
r = '"crit 9ic = 0d- 



4.10 Mass accretion rates 

The rate of stars plunging onto the central SMS is given by two different formulae, depending on whether 
or not there is a crossing point for the 9i c , 9d- plot against the radius, a critical radius. If we find that the 
curves happen to cross, the mass accretion rate has the expression = .-#*(rh)// re iax('"h) because 
the loss-cone will be depleted in a relaxation time, and the mass to take into account is that which lays 
within the critical radius. On the other hand, if there is no crossing point, this means that the loss-cone 
in not empty and in this situation we have to employ a rather different expression, we have to resort 
to the loss-cone star density expression to get the mass being accreted into the SMS. In this case the 
timescale of interest is the dynamical time, Jt, = Q{r)^£*(r) /t& yn {r). 

It may be asserted, nevertheless, that this formula is not completely correct because it is based on the 
stationary model, which supposes an empty loss-cone in the first case and a full loss-cone in the second 
case. We have to generalise it by means of a "diffusion" model (see chapter 5). We introduce the concept 
of the filling degree K of the loss-cone as follows: Let us conjecture that / is the unperturbed velocity 
distribution. If the loss-cone is empty and angular momentum diffusion is neglected, then / = inside 
the loss-cone and / remains unchanged elsewhere in velocity space. Actually, this distribution function 
will have a continuous transition from nearly unperturbed values at large angular momenta towards a 
partially depleted value inside the loss-cone. We approximate this by a distribution function / having a 
sudden jump just at the value L m ; n = m*vi c from an unperturbed value fy, f = K ■ fy, with < K < 1 . 
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Since we work with the hypothesis that some kind of stationary state is to be established in the limit 
t — > oo, the filling degree is the result of calculating this limit, 



(4.44) 



where 



Therefore, 



V(l + V) 



l + v(l + v) 



(4.45) 



(4.46) 



In this expression V = 0^/0^. Then we have to multiply the accretion rates by this filling degree k^. 
The stellar mass within and outside the influence radius is 



M(r,r min )\ r<rh = lim An f p{y J )y jl dr J 

r min^ U •'''mill 

lim 4np(r h )-rl /4 f r> l/4 dr' = 



(4.47) 



lim Anp{r^-rH A 

r min 



5/4 



\6n . , j/rV 

- p(rh)r *u) ' 



M(r)| r>rh -M(r h )| r<rh +4^ /' p(r>' V 



(4.48) 



167T 3 



(l+rVr?) 3/2 (l+^) 3/2 



To get the total heating rates, we just have to compute the value of 



Aieat, all — 1 ^heat, 1*; 



(4.49) 



where £heat, l* is the heating for one star (during one crossing). 
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4.11 Heating rates: An esteem 

In this section we analyse the interaction rate of stars with the SMS by varying the parameters introduced 
in the former sections, namely 8, the core radius r c , the total stellar mass and the super-massive star 
mass itself. We suppose that the stars which conform the stellar system are solar-type stars. For all the 
plots we extend the radii down to 1.001 times the SMS radius, because we would run into a snag if we 
extended it within the SMS radius, since the proportions a 2 oc \/r and p °< r~ 7 / 4 would be wrong and 
the loss-cone star density and therefore the loss-cone angle has been obtained considering a Maxwell- 
Boltzmann distribution. It would also be an error of the problem conception itself, because we are 
studying the non-confined stars and the loss-cone, and its definition does not make any sense for radii 
less than the M. This explains the first inequality in ffl < < r c , which is an exigency that we must 
follow unfailingly because, otherwise, the Plummer model, which we use for the stellar system, would 
not be a suitable solution for it in the case that h, < r c is not satisfied. However, what we demand here is 
not a requirement of the physics of the problem, but a condition for the method being employed to solve 
it. Situations in which rj, < or > r c have to be solved numerically for the equation of Poisson and 
the velocity distributions. 

In order to estimate the heating rates for a single star crossing the SMS we have to plot out the mass 
accretion rate of this central massive object in the case that we have no crossing point for the loss-cone 
and diffusion angle curves, whereas if we have a critical radius we will have to compute the total stellar 
mass and dynamical time at this value. 

In Fig. (4.2) we plot the velocity dispersion against the radius for a 1O 3 M SMS. We observe a typical 
power law of °< r~ x l 2 within the influence radius -which we represent by a vertical dashed line for all 
cases- because we have a cusp on velocities in this interval of radii. We have a nearly constant velocity 
for later radii which lie in the section of values close to the core radius. Then the velocity dispersion 
decays. The length of this nearly constant velocity section as well as the slope for the decay depends on 
the galactic nuclei in the galaxy. 

Regarding Fig. (4.3), we show the loss-cone normalised density difference between the isotropic and 
radially anisotropic cases for a 1O 7 M SMS. In this plot we can observe a bigger number of stars 
being accreted into the SMS for the radially anisotropic case,since a radial orbit means a lower angular 
momentum and thus it is more probable that the star sinks into the central object, and vice-versa for a 
tangential orbit, Fig. (4.4). 

As regards the loss-cone and diffusion angle plots, we examine two different cases in Figs. (4.5) and 
(4.6): a 1O 4 M and a 1O 7 M SMS. One may observe that for the first one we find a critical radius, 
whereas for the latter one the curves do not intersect. It is also interesting to find out which angle is 
bigger and where, since for a 9[ c ^> 0q we have an almost empty loss-cone, because it is quickly depleted 
and the stars replenish it very slowly; the opposite case, 9\ c <C 9d, implies that the loss-cone is full. 
We get a maximum at about a core radius in the mass accretion rates plot for a 1O 7 M SMS, 
because the biggest contribution of stars being accreted into the SMS lies at this radius. Fig. (4.7) shows 
an irregularity at the influence radius, because the loss-cone and diffusion angle plots show us that 
the former happens to be always bigger than this one and thus we have to apply the approximation 
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commented in the foregoing section. 

When we look at the Jl % for the 1O 3 M and 1O 4 M SMSs, we obtain that for the 1O 3 M : ^.| is0 = 
1.75 x 1(T 13 M /yr, J?.\ m = 1.66 x 10~ 13 M /yr, = 1.97 x 10~ 13 M /yr. Where the sub- 

script "iso" stands for isotropy, "rad" for radial anisotropy and "tan" for tangential anisotropy. 
The case of 1O 4 M is similar to the last one. For the 1O 7 M SMS, we will select the corresponding 
to the core radius, since the most important contribution is reached there: ^£,\ com = 10~ 2 M /yr. To 
get the heating rates of these non-confined stars, we just have to compute 



j nrip sms vl ■ 2&/t cmss , (4.50) 

where vl = GJl t cmss = 2&/v+, p SMS = Jt / (^f3$ y ). For the SMS we have supposed, as a first 
approximation, a constant density. 
The corresponding ^"s are: 

^ 10 3 =2.5 x 1(T 9 pc 
J? 10 4 =8x 1(T 8 pc 
M w i =2.5 x l(T 2 pc, 

where stands for the 1O'M SMS radius. 
The luminosities are: 

L lo3 /L =6.2x 10 3 
L W 4/L e = 1.17 x 10 4 
L wl /L & =5.6x 10 s , 

idem L 10 , 1O ! M SMS luminosity. It is not a surprise that these luminosities are not sufficient to support 
quasar luminosities, which was known before. We confirm however the earlier result by Langbein et al. 
(1990) with our more detailed, but stationary loss-cone model, that the luminosities are large enough to 
prevent for some time the relativistic collapse of a SMS in a galactic centre; that was called the quasi-pile 
stage by Hara(1978). 

In the case of a confined star, the heating rate can be obtained in a much simpler way. We just have to 
assume that the star is on a circular orbit at r w then its velocity is 



the orbital time is 



vL = =, (4.51) 



Vci 



arc 



f= 27^> = (452) 



Regarding to the drag force, = Cnnr^p g v C i IC , where p g = ( 4 ^^ 3 ) is the average density. The 
heating by one star will be 

^heating Iconf = 1 2%Sl (4.53) 
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M SMS -10 3 M sun r c = lOpc 5=1 

1 2 1=1 1 1 1 1 1 1 1 1 




2 nil i i I i L i i I i i i i 

-5 5 

Log(r) 



Fig. 4.2: The velocity units are (m/s) 2 and r is expressed in pc. The influence radius is located at 6- 10~ 4 pc. 
We show its logarithm in the plot with a vertical dashed line. The central velocity dispersion (the velocity 
dispersion at the influence radius) is 0" centra i = y/ GJ?* / {6r c ) = 84 km/s. Msms stands for J( 



The number of stars in the SMSis given by 

N*=^-p(0)&, (4.54) 
where p(M) = p(t\)(M Ah) - ' ?/ ' 4 '■ Therewith, the heating by confined stars is 

^heating, conf = N* ■ F D ■ 2%3%. (4.55) 
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M SMS =10 7 M sun r c = lOOpc 5 = 10 




Log(r/r c ) 



Fig. 4.3: Q(r) := Pic( r )/P( r ) decays with the distance to the SMS. We have divided the radius by the core 
radius in order to normalise it. A bigger number of stars plunge on to the SMS for the radially anisotropic 
case 

4.12 Discussion of the results 

We have revisited the classical loss-cone semi-analytic theory invented by Frank and Rees (1976) for 
star accretion onto central super-massive black holes in galactic nuclei and star clusters and extended by 
da Costa (1981) for the case of stars on radial orbits being trapped by star-gas interactions in a central 
super-massive star-gas system. In Langbein et al. (1990) such model was included in time-dependent, 
spherically symmetric models of star-gas systems in galactic nuclei. Though highly idealized, we think 
that such configurations are still worth a study to understand the physical processes at work in the early 
formation phase of massive galaxies with formation of central black holes. Notions such as the critical 
radius in the classical work, where the loss-cone star accretion becomes important and flattens out the 
cusp density profile turn out naturally in our model without any ad hoc assumptions. While in this 
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M SMS =10 7 M sun r c = lOOpc 6 = 10 




-2-10 1 2 

Log(r/r c ) 



Fig. 4.4: The same as for Fig. (4.3) but for the tangentially anisotropic case and a lower number of stars 



research note we keep the stellar background system fixed and develop the old ideas in an up-to-date 
form, one should include them into a self-consistent dynamical model of relaxing star clusters with 
central black hole. Recent dramatic improvements of the observational situation and the demography 
of black holes in nuclei demand such progress in dynamical modelling, which is yet surprisingly poor 
(Richstone et al., 1998). It is clear that in a proper cosmological context many complications occur 
for simplified modelling: the whole system is embedded in a collisionless dark halo with disputed 
central density profile, it is non-stationary due to a sequence of merger events in hierarchical structure 
formation, with the possible formation of binary or multiple black holes and perturbations of various 
kinds will even cause a single black hole not to be fixed in the centre. Despite of all that have begun 
our work at the point where we think present (astro)physical understanding and modelling comes to its 
limits, and this is the case for a spherical dense large N star cluster, suffering from relaxation and star 
accretion around a fixed massive black hole. 
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Fig. 4.5: &i c =2 - yj n(r), 6*d = \AdynArelax- For this mass the two curves cross at the critical radius. For 
radii smaller than that, the loss-cone angle is bigger than the diffusion angle, and this implies that the 
loss-cone is empty. From the r cr j t onwards it is no longer empty, for > 0[ c . With the crossing point we 
can work out the accretion rate and find out the differences depending on whether we consider an isotropic 
situation for the stellar velocity distribution function or an anisotropic one, distinguishing between a radial 
or tangential anisotropy. If we set this case against the 10 3 M e we do not find big differences 



4.13 Addendum A: The tidal radius 

When the stars approach the BH they experience strong tidal forces which may be violent enough to 
blow up the star. The star gets disrupted whenever the work exerted over the star by the tidal force 
exceeds its own binding energy, (all energies are per unit mass): 

£bind = a , a = (4.56) 
r± 5-n 

where n is the polytropic index (Chandrasekhar, 1942). 

( Fl -F 2 )2r* = a^ (4.57) 
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M SMS =10 7 M sun r c = lOOpc 5=10 




-3-2-10 1 2 

Log(r/r c ) 

Fig. 4.6: For higher masses, such as 10 7 M S , the loss-cone is empty. We get no critical radius 



GJZ. p _ GJZ, 

(r t -r*) 2 ' (n + r+) 

Considering r* <C r t , we can approximate the expressions: 



^ = 7—T-2^2^ T —- V2 . (4.58) 



1 , 2r, 



(r,-r,)2 



&h (4 - 59) 



then, 
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M SMS =10 7 M sun r c = lOOpc 5=10 




-2 2 

Log[r/r c ] 



Fig. 4.7: The "jump" is due to the approximation we made in order to get the accretion rates. The maximum 
is reached at a distance which corresponds to the core radius 



'"tidal 



-(5-n) 

3 



1/3 



(4.60) 



For solar-type stars it is (considering a n = 3 polytrope) 



r t ~1.4x 10 1 



Mr. 



1/3 



cm. 



(4.61) 
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Fig. 4.8: Decomposition of the tidal forces over a star 



4.14 Addendum B: The difussion and loss-cone angles 



4.14.1 Definition of the difussion angle 9d 

The relaxation time is the required time for Av^/v^ ~ 1 (i.e. the change in the perpendicular velocity 
component is of the same order as the perpendicular velocity component itself); 

Avi = » relax • Svi, Avi/vi = 1 = " relax 2 gv ' (4.62) 

( v\ \ 

^relax — n relax ' fdyn — | ) ' ^dyn? (4.63) 



where n re i ax is the numbers of crossings for A v\ / v]_ ~ 1 . This conforms the definition of the relaxation 
time, Av^/v^ = f/f re iax (Binney and Tremaine, 1987). 
If we consider that D is very small, 



— 0D> 'relax — 2 

v y n 



sin 9 D ~ ~ D ; f re iax - (4.64) 



9d^/-^. (4.65) 

'relax 
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V 




Fig. 4.9: Definition of angle 




• BH 

Fig. 4.10: Definition of angle 6\ c 

4.14.2 Definition of the loss-cone angle 0\ c 

We suppose that the central object with mass Ji % has an influence radius Hi. To define this radius we 
say that a star will interact with the central object only when r < i\. Then we look for a condition at 
a place r > for a star to touch or to cross the influence radius of the central object within a crossing 
time t cross = r/a T . 

As we saw in the text, the condition that defines this angle is the following: 

r p (E,L) <r u 6< 6*. (4.66) 

Vt Vt L/r 

sin0 = -, 0<1; =*0~- = -i-. (4.67) 

V V V 
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Where L — rv t is the specific angular momentum. 

Now we derive an expression for this angle in terms of the influence radius: 
• r < Hi 

In this region the star moves under the BH potential influence, then 



a (r) „ ^ 9^L = ^±^. ^ = a{Rh) . ^f r = ac . v /^77, (4.68) 
since a t 2 = G^#.//?h- 

The typical velocity of the orbit is < v 2 >~ 3 a 3 , where the factor three stands for the three directions in 
the space. Since a means the one- dimensional dispersion, we have to take into account the dispersion 
of the velocity in each direction. Then, 

<v>~^/3G c ^/nJ~r. (4.69) 

Finally, we obtain the loss-cone angle, 

V 3 r 

• r>r b 

We can consider that the velocity dispersion is more or less constant from this onwards, v w V3o c , 



6 k = V2 ^' rx ;o c = ^G^. (4.71) 



The angle is 



0k- ~ \l\r trh /r. (4.72) 



4.15 Addendum C: The loss-cone velocity 

We derive the v\ c {r) using angular momentum and energy conservation arguments. We just have to 
evaluate it at a general radius r and at the tidal radius r u where the tangential velocity is maximal and 
the radial velocity cancels (see Fig.4.1 1). 
• general radius: 

L(r)= rv tg {r) (4.73) 
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Fig. 4.11: Definition of the peribarathron as the distance of closest approximation of the star in its orbit to 
the BH. In this point the radial component of the velocity of the star cancels and the tangential component 
is maximum. In the figure "r p " stands for the peribarathron radius and "r t " for the tidal radius 



• tidal radius: 



L{r t )= r t v tg (r t ), (4.74) 



from the momentum conservation and the fact that v r (r t ) = 0we get that: 



vtg(r t ) = f -vtg(r). (4.75) 



Idem energy, and using the last result, 



2 „ (r\2 r 2 



Vtg(rY v r (r) 



Hr) - ~ = Hn) ~ ^ ■ v tg {rf. (4.76) 

Then we get the tangential velocity of the stars in terms of r; namely, the loss-cone velocity: 
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Therewith, the angular momentum is 



Vk M = -1^=2 ■ \/2Mn)-<l>(r)]+Mr) 2 . (4.77) 



r 

L(r t )= r t -v t g(r)| max = r,--vtg(r) 

ft 



r 



where 



If we use the fact that 



v tg (r) =r Zl_^ ^2^ + v v {r)\ (4.78) 



A</> = </)(r t ) - 0(r) = + ^(r t ) - ^(r) (4.79) 



r>^_»(_ + *(r)j=,(r) (4.80) 

Also it is 

> 0*(r t ), since > ^f*(r t ) (4.81) 

r t 

thus, 



r V n 



If we use now the fact that 



^_» hH i)- w _^.(i)- ,fl 



'^.^.flV" 1 (4,4, 



It is then in fair approximation 



vic(r)«-^4/-jp«^(r).M] . (4.85) 
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A diffusion model for accretion of stars 



IEWING the problem of the joint evolution of a spherical star cluster with a central BH making 



feasible anisotropy, we shall introduce in these pages the simulation method to follow to anal- 



T yse such a scenario. This procedure, as we exposed in former chapters, is based on moments 
of the Boltzmann equation with relaxation. The cluster is modelled like a self-gravitating, conducting 
gas sphere, according to the methods presented in Louis and Spurzem (1991) and Giersz and Spurzem 
(1994). These models improved earlier gas models of Bettwieser (1983) and Heggie (1984). Much as 
the structure of the numerical method is for the sake of computational efficiency on account of physi- 
cal accuracy, it allows for all the most important physical ingredients that may carry out a role in the 
evolution of a spherical cluster. These include, among others, self-gravity, two-body relaxation, stellar 
evolution, collisions, binary stars etc and, undoubtedly, the interaction with a central BH and the role of 
a mass spectrum. The specific advantage of the so-called "gaseous model" to other simulations methods 
(see chapter 2 for a description) is that the simulations are comparatively much faster, since they are 
grounded on numerical integration of a relatively small set of partial differential equations with just one 
spatial variable, the radius r. In addition, all quantities of interest are accessible as smooth functions of 
r and time and this allows one to investigate in detail clean-cut aspects of the dynamics without being 
hindered by the important numerical noise particle-based methods (N-body and Monte Carlo) suffer 
from. 

In this chapter we concentrate on the simplest version of the gaseous model which includes the interac- 
tion between a central BH and its host cluster. In particular, we assume that all stars are sun-like, neglect 
stellar evolution, direct collisions between stars and the role of binary stars. Also, the only interaction 
between BH and the stellar system is tidal disruptions (besides the BH's contribution to the gravita- 
tional field), and we undertake that the BH stays fixed at the centre of the cluster. While admittedly 

1 We employed some of the results of this chapter for Amaro-Seoane et al. (2004) and Spurzem et al. (2003) 
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very simplified, we reckon this idealised situation warrants consideration. First, it helps us to establish 
that the gaseous model is also able to treat more complex situations in phase space than self-gravitating 
star clusters, such as those caused by loss-cone accretion and a central BH; second, experience with the 
gaseous model and other methods showed us how intricate the interplay between various physical ef- 
fects can become during the evolution of clusters and so we feel compelled to first consider the simplest 
models to develop a robust understanding of the mechanisms at play. 

The structure of this chapter is as follows: In section 5.1 we explain the physics of the problem. In 2.2 
we introduce briefly the various analytical and numerical methods used so far to investigate spherical 
clusters with central BH and summarise their key results. In section 5.2 we give a description of the 
theoretical diffusion model. Section 5.4 is devoted to results obtained in a first set of simulations. 
Finally, in 5.5, we draw conclusions about this first use of this new code version and present what our 
future work with it will likely consist of. 

5.1 Loss-cone accretion on to massive BHs 
5.1.1 Previous theoretical and numerical studies 

If we arrange numerical methods for stellar dynamics in order of validity and both increasing the spatial 
resolution and decreasing the required computational time, we can distinguish four general classes. The 
most direct approach is the so-called A^-body method (Aarseth, 1999a,b; Spurzem, 1999). Monte Carlo 
codes are also particle-based, but rely on the assumptions that the system is spherically symmetric and 
in dynamical equilibrium and treat the relaxation in the Fokker-Planck approximation (see 3.2) (Freitag 
and Benz 2004; Fregeau et al. 2003; Freitag and Benz 2001; Giersz and Spurzem 2003; Giersz 1998; 
Joshi et al. 2001). Ensuingly, we have the two-dimensional numerical direct solutions of the Fokker- 
Planck equation (Takahashi, 1997, 1996, 1995), and the gaseous models. The idea of this model goes 
back to Hachisu et al. (1978) and Lynden-Bell and Eggleton (1980), who first proposed to treat the two- 
body relaxation as a transport process like in a conducting plasma. They had been developed further 
by Bettwieser (1983); Bettwieser and Sugimoto (1984); Heggie (1984); Heggie and Ramamani (1989). 
Their present form, published in Louis and Spurzem (1991); Giersz and Spurzem (1994); Spurzem and 
Takahashi (1995) improves the detailed form of the conductivities in order to yield high accuracy (for 
comparison with iV-body) and correct multi-mass models. This point has been made already in Spurzem 
(1992). 

Peebles (1972); Shapiro and Lightman (1976) and especially Frank and Rees (1976) and Bahcall and 
Wolf (1976) addressed the problem of a stationary stellar density profile around a massive star accreting 
BH. They found that, under certain conditions, the density profile p °< r~ 7 / 4 is established in the region 
where the BH's gravitational potential well dominates the self-gravity of the stars 2 . 
The problem of a star cluster with a massive central star-accreting BH has been widely coped with 
Fokker-Planck numerical models. This approach was useful in order to test the solidness of the method 

2 We must mention here the legwork done twelve years before this analysis by Gurevich (1964), since he got an analogous 
solution for the distribution of electrons in the vicinity of a positively charged Coulomb centre. 
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to reproduce the p r~ 7 / 4 stationary density profile, since loss-cone accretion disturbs such a density 
cusp (Ozernoi and Reinhardt, 1978). The authors show that the stationary density profile follows from 
their stellar-dynamical equation of heat transfer by scaling arguments which are analogous to those 
given in Shapiro and Lightman (1976). 



5.1.2 The diffusion model 

We can express the tidal radius in terms of the internal stellar structure re-writing Eq. (4.60), 



X b 7tp 



(5.1) 



where denotes the mass of the central BH, p the mean stellar internal density, n is the polytropic 
index (stars are supposed to be polytropes) and x\, is a parameter proportional to the gravitational binding 
energy of the star that describes effects of the internal stellar structure. We assume that a star is disrupted 
by tidal forces when it crosses the tidal radius. The free parameter e e ff (accretion efficiency) determines 
the mass fraction of the gaseous debris being accreted on to the central BH (e e ff = 1 corresponds to 
100% efficiency). 

There are two concurrent processes driving stars towards the tidal radius; namely the energy diffusion 
and the loss-cone accretion. In the first case, stars on nearly circular orbits lose energy by distant 
gravitational encounters with other stars and in the process their orbits get closer and closer to the 
central BH. The associated energy diffusion time-scale can be identified with the local stellar-dynamical 
relaxation time, Eq. (1.1), but generalised for anisotropy as in Bettwieser (1983), 

reax \6y/n G 2 m- k p+{r) In A 

Here c r , a t are the radial and tangential velocity dispersions (in case of isotropy 2a r 2 = of), p*(r) is the 
mean stellar mass density, N the total particle number, G the gravitational constant, m+ the individual 
stellar mass and 



In A = In Cp max //?o) = In (yN) (5.3) 

is the Coulomb logarithm. We set y = 0.11 (Giersz and Heggie, 1994). In this expression p max is an 
upper limit of p, the impact parameter, po is the value of p that corresponds to an encounter of angle 
y/ = n/4, where y/ = (n — §)/2 if £ is defined to be the deflection angle of the encounter (Spitzer, 
1987). In the vicinity of the BH (r < r^, see Sec. 5.3), one should be aware that A w ^t,jm* (Bahcall 
and Wolf, 1976; Lightman and Shapiro, 1977) but, for simplification, here we shall use Eq. 5.3 (strictely 
speaking only valid at distances r > Hi) everywhere. 

For a more detailed discussion of the energy diffusion process and its description in the context of the 
moment model see Bettwieser and Spurzem (1986). 

As regards the second process, the loss-cone accretion, stars moving on radially elongated orbits are 
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destroyed by tidal forces when they enter the tidal radius r t . A star will belong to the loss-cone when 
its peribarathron (distance of closest approach to the BH, see Fig. 4.1 1) is less than or equal to the tidal 
radius r t , provided that its orbit is not disturbed by encounters. Thus, the loss-cone can be defined as 
that part of stellar velocity space at radius r, which is given by 

\v t \< Vfc(r) = -^jL= ■ J2{<p(r t ) - <p(r)]+v r (r) 2 (5.4) 

V r - r t 

(see Eq. 4.77). In the last formula v r , v, are the radial and tangential velocity of a star and 



0(r t )-*(r) = - + Wt)- 




(5.5) 



At distances r>r t we can approximate this expression taking into account that G^,/r t 3> (r) and (r t ) 
(seeEq. 4.85). 

For a deeper analysis on loss-cone phenomena see chapter 4. 

Similar to da Costa (1981), we define two time-scales: f out to account for the depletion of the loss-cone 
and f; n for its replenishment. For a BH, f out is equivalent to one crossing time, since it is assumed that a 
star is destroyed by just one crossing of the tidal radius. In chapter 4 we mentioned that in the special 
case that the central object is a super-massive star, f out = nt cmss . Here n > 1 is the number of passages 
until a star is trapped in the central object. 

The loss-cone is replenished by distant gravitational encounters that change the angular momentum 
vectors of the stars. To estimate f; n , we make the assumption that random gravitational encounters 
thermalise the whole velocity space at some given radius r after a time-scale of the order of f re iax- As 
a first approximation, the fraction Q. of the three-dimensional velocity space is re-populated within a 
time-scale 



f; n — Q. f re lax (5.6) 

where 

Q:= / fd 3 v/ [ fd 3 v; (5.7) 

the subscript °° denotes an integration over the velocity space as a whole, and the subscript Ic means 
that the integration over the loss-cone part of the velocity space is given by Eq. (5.4). As a matter of 
fact, Q. can be envisaged as the fraction of the surface of a velocity ellipsoid which is cut out by the loss- 
cone. Close to the tidal radius r t , and for appreciable amounts of stellar-dynamical velocity dispersion 
anisotropics, our method describes the loss-cone size Q. more exactly than those given by (Frank and 
Rees, 1976; da Costa, 1981). On the other hand, their models can be "recovered" in the limit of r ^> r t 
and isotropy, 2 of w of (denoted as the "small loss-cone approximation"), where 

i2«vf c /o t 2 «r t /r. (5.8) 
This is equivalent to their definition of d\ c if Q. = 0^/4 is adopted. 
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Where the loss-cone effects can be neglected, a Schwarzschild-Boltzmann type distribution function can 
be assumed, 

Third order moments of the velocity distribution that represent the stellar-dynamical energy flux do not 
alter such a distribution function significantly (Bettwieser and Sugimoto, 1985). 

An important quantity is the critical radius r cr j t . Let 9^ be the average quadratic deflection angle 
produced by relaxation during r out (= f cr0 ss here), 

0g = i2HL. (5.10) 

frelax 

Then, by definition, 

which is equivalent to f ut(r c rit) = 4/i n (r cr it). For most clusters where such a radius can be defined, 
01c > 0d inside r cr it while the opposite holds outside. This means that, at large radii, relaxation is 
efficient enough to make stars diffuse into and out of loss cone orbits over a time scale f ou t so that the 
distribution function is not appreciably depleted in the loss cone. Conversely, deep inside r cr i t , the loss 
cone orbits are essentially empty and the flux of stars into this domain of phase-space (and into the BH) 
can be treated as a diffusive process because the size of one individual step of the velocity random walk 
process, 0q, is (much) smaller that the characteristic size of the problem, 9\ c . 

Note that a critical radius does not necessarily exist (see chapter 4). For instance, if one assumes that 
gravity of the BH dominates the stellar self-gravitation and that the density profile follows a power-law, 
p cx r~ a , one has 0^ « 0^ °< r 3 ~" and a critical radius would not exist for a > 4. 
Now we want to generalise the stationary model (see chapter4), which assumes an empty loss-cone 
within r cr jt and a full loss-cone elsewhere, by means of a simple "diffusion" model, which is derived 
from the above considerations; this means that the filling degree of the loss-cone K can be continuously 
estimated within its limiting values, 

Ke [0,1]. (5.12) 

Let / be the unperturbed velocity distribution (without loss-cone accretion); if the loss-cone is empty 
and we neglect the angular momentum diffusion, / = inside the loss-cone (and unchanged elsewhere in 
the velocity space). In point of fact, / will have a continuous transition from nearly unperturbed values 
at large angular momenta to a partially depleted value within the loss-cone. This value is determined 
by the ratio of f; n and f out . Such a smooth transition of the distribution function given as a function of 
angular momentum /(/) has been derived from self-consistent models of angular momentum diffusion 
(e.g. Cohn and Kulsrud 1978 or Marchant and Shapiro 1980). We approximate f(J) by a distribution 
function that has a sudden jump just at the value J m i n = m*v\ c from an unperturbed value fo given by the 
moment equations (assuming a Schwarzschild-Boltzmann distribution) to the constant lowered value 
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f = Kf , withO<ZT< 1 (5.13) 

within the loss-cone (i.e. J < 7 m ; n or \v t \ < vi c ). This implies that, as means to compute the mean mass 
density of loss-cone stars, we have to calculate the integral 

pi c = f Kf d\. (5.14) 

Jlc 

And then, accordant with the definition of il, 

Pic,faU=P^, (5.15) 

in the case that we have a full loss-cone. 

In regard to the radial and tangential stellar velocity dispersions in the loss-cone a\ c r and Oi c r , we can 
compute them using second moments integrated over the loss-cone part of velocity space. As for the 
definition of the quantities E r and E, used in Sect. 3, 

°ic,r — £r°r i 

<, - E t ol (5.16) 
in the small loss-cone approximation we have that £ r wl and £,<1. 

The arguments about the time-scales that have led us to the derivation of t m and f ut guide us also to the 
following diffusion equation for the time evolution of the spatial density pj c = Kp Q. of loss-cone stars: 

dp\c = picflc | pfl-Pk (5 1?) 

dt f ou t t{ n 

In this equation, the second term on the right hand is the refilling term due to relaxation. 
As we assume relaxation is due to a large number of small-angle deflections and can thus be seen as a 
diffusive process in velocity space, the probability P(9) that a star is scattered in an angle 9 in a time 
t ov x is 



P(0) = — — exp(-0 2 /02 ); (51g) 



The distribution is normalised to one, 







P(6)d6 = \ (5.19) 



and has the property that its mean square value is 0^. A star remains in the loss-cone during a time ? out 
if its RMS diffusion angle is smaller than the loss-cone angle 9i c . The probability for this to happen is 

P lc = / P(0) C /0=erf(v^n7w). (5.20) 
Jo 
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In the case that a star is unperturbed by the rest of the stellar system, it will sink on to the central BH 
in a time f ou t- Actually, this is a somehow simplified description of the physical process, for part of the 
loss-cone stars will be scattered out of it before they slump. The required time for this event is f; n , since 
in this time-scale the angular momentum vector will change (due to distant encounters) on an amount 
that is comparable with the size of the loss-cone in the angular momentum space. For this reason we 
have introduced the quantity F\ c in Eq. (5.17). 

Bluntly speaking, the effective time-scale that describes the loss-cone depletion allowing for perturbation 
due to angular momentum diffusion is 

fout.eff = fout/flc- (5-21) 

As a matter of fact, this definition ensures us that far outside of the critical radius the loss-cone depletes 
in a time that grows infinitely, as it is physically expected. In the regime where r <C r c rit, P\ c tends 
asymptotically to 1 and to where r 3> r c rit> passing through a transition zone at r = r cr j t . 
We can consider Eq. (5.17) as an ordinary differential equation for K = p\ c / (p Q.) if we assume that the 
stellar density and the loss-cone size are time-independent. Transport phenomena can be neglected, for 
they are related to the relaxation time, and f re i ax 3> fj n , f out . Bearing this in mind we can get an analytical 
solution K(t) for the differential equation with the initial condition that K(t) \ tf) = Kq, 



v 'out ' 

^•(— (-^)> 

In the last equation we have defined E, for legibility reasons as follows, 



4 := l + (Wflc«in). (5.23) 
For r = r cr it, with f; n = r out , the stationary filling degree of the loss-cone turns out to be 

K oa :=]imK(t) = l. (5.24) 

Note that Milosavljevic and Merritt (2003) recently gave a detailed summary of loss-cone effects. They 
derived expressions for non-equilibrium configurations. They employ a rather different treatment for 
the diffusion since they tackle the problem of binary BHs scattering. 



5.2 Inclusion of the central BH in the system 

In this subsection we discuss the way we cope with the loss-cone in our approach. For this aim we 
accept the following: 

1 . The system has central (r = 0) fixed BH 
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2. Stars are totally destroyed when they enter r t 

3. Gas is completely and immediately accreted on to the BH 

As regards the first point, one should mention that the role of brownian motion of the central BH can 
be important; as a matter of fact, for a cluster with core radius R CO k, equipartition predicts a wandering 
radius /? W an of order 



flwan W flcore y/M*/^., (5.25) 

which is larger than the tidal disruption radius for BHs less massive than 1O 9 M if the core radius is 
1 pc (Bahcall and Wolf, 1976; Lin and Tremaine, 1980; Chatterjee et al., 2002). The wandering of a 
MBH at the centre of a cuspy cluster has been simulated by Dorband et al. (2003) with a A^-body code 
allowing N = 10 6 . They find that RMS velocity of the quickly reaches equipartition with the stars but 
do not comment on the wandering radius. In Addendum A, at the end of this chapter (section 5.6), 
we present a simple estimate suggesting that, in a cusp p oc r ~ a , the wandering radius may be much 
reduced, /? wan °< a(m+/ ^.) l K 2 ~ a \ where a is typical length scale for the central parts of the cluster. 
For a > 1.5, one would then expect R wan to be smaller than R t for black holes as light as 2000 M , 
but this arguments neglects the flattening of the density profile due to loss-cone accretion. Further N- 
body simulations are clearly required to settle the question and, in particular, if /? wan is larger than R t , 
to establish the effect of the motion of the MBH on disruption rates, which can be either increased or 
decreased (Magorrian and Tremaine, 1999). 

We arrogate that at any sphere of radius r the transport of loss-cone stars in the time-scale f ut,eff towards 
the centre happens instantaneously compared with the time step used for the time evolution. 
Hence, the local density "loss" at r is 



<Sp\ Pi A , 

<St/ lc tout 



(5.26) 



with pic = KQ.p. 

This corresponds to a local energy loss 3 of 



8t Ac \ St Jlc 



E t al (5.27) 



E r and E t are worked out integrating over the velocity distribution part that corresponds to the loss-cone 
with the approximation u <C a r and vi c <C C7 t , 



3 By loss we mean here transport of mass and kinetic energy toward the central BH, for it is lost for the stellar system. 
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E r K, 1 

E t w vf c /cJ t 2 <l (5.28) 



Thereupon, the mass accretion rate of the central BH can be calculated as 



8p 



M=-e eB J 4nr 2 dr. (5.29) 

Here Rm stands for the total radius of the stellar system. The accretion efficiency has been set throughout 
our calculations to £ e g = 1; for a discussion on different e e ff-values see Marchant and Shapiro (1980). 
The complete set of Eqs. (2.22) including the local accretion terms of the type (8/8t)i c for energy 
diffusion and loss-cone accretion are solved implicitly. For every time step the mass of the BH and the 
filling degree K of the loss-cone are brought up to the new state of the system. The time step is chosen 
in order to keep the maximum changes of the variables below 5% . 

For the model calculation we have utilised for the boundary conditions that at the outer limit, R tot = 
10 4 pc, we impose u = F = and M r = M tot . No stellar evaporation is allowed. At the centre, the 
usual boundary conditions for the gaseous model are u = M v = F = but the central point r\ = is 
not explicitly used when there is a BH, for obvious reasons. Instead, one imposes that all quantities 
vary as power-laws, d lnx/d In r = C st inside the first non-zero radius of the discretisation mesh, r 2 = 
1.7 x 10~ 6 pc (see chapter 2). 



5.3 Units and useful quantities 

The units used in the computations correspond to the so-called N-body unit system, in which G = 1, 
the total initial mass of the stellar cluster is 1 and its initial total energy is —1/2 (Henon, 1971; Heggie 
and Mathieu, 1986). For the simulations presented here, the initial cluster structure corresponds to 

the Plummer model whose density profile is p(r) — po( 1 + (r/Rpi) j , where R P i is the Plummer 
scaling length. For such a model the A^-body length unit is = 16/(3?r)/?pi. 

In the situations considered here, the evolution of the cluster is driven by 2-body relaxation. Therefore, 
a natural time scale is the (initial) half-mass relaxation time. We use the definition of Spitzer (1987), 

For a Plummer model, the half-mass radius is R^ 2 = 0.769"^ = 1.305/?pi. *^ci is the total stellar mass. 
For a cluster containing a central BH, an important quantity is the influence radius, enclosing the central 
region inside of which the gravitational influence of the BH dominates over the self-gravity of the stellar 
cluster. The usual definition is r^ = GJi % j C7q , where Gq is the velocity dispersion in the cluster at a large 
distance from the BH. As the latter quantity is only well defined for a cluster with an extended core, we 
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use here the alternate and approximate definition M r (rh) = i.e. rj, is the radius of that encloses a 
total stellar mass equal to the mass of the BH. 

5.4 Results 

We study the evolution of a stellar cluster with a so-called "seed BH" at its centre. We consider two 
possible configurations for the stellar system; one of a total mass of M tot — 1O 5 M and another of 
1O 6 M . For the initial BH mass, we have chosen ^#.(0) = 50M Q and 500M© and we model it as a 
Plummer of R P \ = lpc. Even though it would be more realistic to set the efficiency parameter e e s =1/2, 
we choose here e e ff = 1 for historical reasons. Nonetheless, here we study additionally the influence of a 
variation of the stellar structure parameter X},, since it influences the tidal radius and hence the accretion 
rates (see Eq. 4.60). For this intent we compare case a (xi, = 1) with another one in which we choose the 
value Xb — 2 (case b). As regards the physical meaning, the stars of case b have twice as much internal 
binding energy than case a. 

The cluster evolves during its pre-collapse phase up to a maximum central density from which the 
energy input due to star accretion near the tidal radius becomes sufficient in order to halt and reverse 
the core collapse. Immediately afterwards, the post-collapse evolution starts. At the beginning of the 
re-expansion phase, the BH significantly grows to several 10 3 solar masses. Thereon, a slow further 
expansion and growth of the BH follow. 

In Fig. 5.1, we follow the evolution of the mass of a central BH in a globular cluster of 10 s stars of 
1 M . Panel (a) shows the mass of the BH as function of time. On panel (b), we present the accretion 
rate on to the BH, i.e. its growth rate. For .^bh(0) = 50 M , the early cluster's evolution is unaffected 
by the presence of the BH which starts growing suddenly at the moment of deep core collapse, around 
T ~ 14.5 7i-h(0). In Fig. 5.2 we follow the same evolution for the case of a stellar cluster of 10 6 stars. 
From Figs. (5.1) and (5.2) we can see that the differences between the cases a, b and c are nearly 
negligible after core collapse. In general, the structure of the cluster at late times is nearly independent 
of ^.(0) and xt,. From these plots we can infer that this occurs since core collapse leads to higher 
densities if the initial BH mass is smaller and thus the integrated accreted stellar mass increases. 
We exhibit the evolution of the structure of the cluster for case a with 10 stars in Fig. 5.3. With dotted 
lines we plot various Lagrangian radii, for mass fractions ranging between 10~ 3 % (which formally 
corresponds to only one star) to 90%. Only the mass still in the stellar component at a given time is 
taken into account. Moreover, the evolution of the influence radius (solid line, defined as the radius 
enclosing a stellar mass equal to the BH mass) and critical radius (dashed line) are shown, so that one 
can infer the percentage of the stellar mass embodied within them at a certain moment. For late time, 
one obtains self-similar evolution with size increasing like R °< T 2 / 3 , as expected for a system in which 
the central object has a small mass and the energy production is confined to a small central volume. 
(Henon, 1965; Shapiro, 1977; McMillan et al., 1981; Goodman, 1984). We consider too the case of a 
10 6 stars in Fig. 5.4, for which the R °c r 2 / 3 expansion is a poor approximation because, at late times, 
the BH comprises of order 40 % of the system mass. 

We observe in Fig. 5.5 that for the evolved post-collapse model the spatial profile of the stellar density 
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Fig. 5.1: Evolution of the mass of a central BH in a globular cluster of 10 5 stars of 1M©. We considered 
three cases. In case a (solid line), the initial BH mass is ^.(0) =50M© and x/, = 1, case b (dashes) has the 
same initial BH mass but jq, = 2 while case c (dash-dot) corresponds to ^#.(0) = 5OOM and x/, = 1. An 
accretion efficiency of e e ff = 1 is assumed. Panel (a) shows the mass of the BH as function of time and panel 
(b) the accretion rate on to the BH. At late times, the mass of the central BH increases like jM m °= r -1 - 2 as 
predicted by simple scaling arguments (see text) 
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Fig. 5.3: Evolution of the radii of spheres enclosing the indicated fraction of the total mass, Lagrangian 
radii, for case a. The mass fractions range from 10~ 3 % to 90%. The influence and critical radius are 
displayed (solid and broken line). See text for further explanation. 



has a power law slope of p « r~ 7 / 4 in the region r cr j t < r < r\, where i\ is the influence radius. The 
density profile flattens for r < r cr j t due to the effective loss-cone accretion. For the same post-collapse 
moment we display the surface density for case a in Fig. 5.8. 

Arguments based on constant accretion rate lead to the result that one should expect a stellar density 
profile proportional to r~'/ 2 inside r cr i t . If one takes Eq. (5.17) and considers that M = C st , we get that 

%ocr-«- 5 / 2 . (5.31) 
at 

Where we have assumed that F\ c as 1, p °c r~ a and suppose that the gravity of the BH totally dominates 
this region, so that we have a Keplerian profile, a °^ r~ x l 2 . We have also taken into account the fact that 
fdyn = r/a. Thus, 

M = r 1/2 - a . (5.32) 
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Whence, a = -1/2. 

Another possible way to derive this result is based on the Eq. (4-137) of Binney and Tremaine (1987), 

p (r) = 4tt [ ¥ f(e)^2( W ~e)de (5.33) 
Jo 

In this equation y := — + 0o and e := — E + <po, where <po is a constant that makes / > for e > and 
/ = for £ < 0. If we assume that f (e)=0 for £ > £o (where £o is around r cr ; t ) and consider that we are 
in the region where r <C r cr i t , y/ = G^£,/r ^> £ and then, from 

p{r) xiAnJly / f(e)de (5.34) 
Jo 

one can conclude that p{r) °^ r~ l l 2 . 

Nevertheless, none of these derivations can be applied to the interpretation of our results, since it would 
contradict the basis of our model's idea. We cannot expect this "classical" result in our results, since it 
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Fig. 5.5: Density profile for a 10 5 M Q globular cluster with ^,h(0) = 5OM (case a) at 10 Gyrs. The round 
dot indicates the influence radius (M r = ./#bh). the star the critical radius and the square the radius below 
which the description of the cluster as a continuum loses significance because the enclosed mass is smaller 
than lMg. As expected, for the zone between the critical and influence radii, the density profile closely 
reassembles a power-law of exponent —7/4. At that stage, the structure of case c (^bh(O) = 5OOM ) is 
extremely similar. We can see that from the "1-star" radius onwards the slope of the curve shows a tendence 
to increase. This is due to the fact that the loss-cone size is artificially limited for stability purposes. On 
the other hand, the slope comprised between the critical radius and the 1-star radius is consistent with the 
arguments given in the work of Lightman and Shapiro (1977) (see text for further explanation). 



would imply that the filling degree K is constant in the region where fl c «l(r« r cr i t see Fig. 5.12), 
which is not the case at all. 

As a matter of fact, one does not know a priori which kind of slope one will have in this region. It is 
rather an initial condition that we introduce and it is totally disconnected with our diffusion model. We 
could bargain for such a slope in the case that we just have stars that come from unbound orbits and fall 
on to the central BH disappearing immediately from the stellar system. In the model presented here, we 
have also stars with bounded orbits around the BH that will perturbe the slope of a = —1/2. Lightman 
and Shapiro (1977) proved that within r cr j t , according to their Eq. (71) (where they assume £2 <C 1, small 
loss-cone), a will continously vary from 1.75 to — °o. 

Regarding the three dimensional velocity dispersion, in Fig. 5.6 we can see coming a slope of -1/2 at 
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T= 9.63x 1 0' F rh (0) = 1 .00x 1 10 yrs 
Radius in pc 
0.1 1 10 10 2 




Radius in N-body units 

Fig. 5.6: Profile of the three-dimensional velocity dispersion for the same case and time as Fig. 5.5. See 
text for comments. 



the inner region; but extending below the one-star radius does not make much sense. A slope of -1/2 
is what one would expect from Kepler's third law and a simple application of Jeans equation, with 
the assumptions that (1) dynamical equilibrium holds, (2) the gravity is dominated by the central BH, 
(3) the density follows a pure power-law and (4) the anisotropy a t /<J r , is constant, indeed predicts 
(7 oc r~ x l 2 . At the end of this chapter, in Addendum B (section 5.7) , we show that the Jeans equation for 
stationary equilibrium actually describes the central regions of the cluster quite well. The reason why 
the velocity dispersion does not follow closely the "Keplerian" profile has to do with the fact that none 
of assumptions (2)-(4) exactly holds all the way from the influence radius inward. 
Figures 5.9 and 5.8 give the plots of the projected density and velocity dispersions for the late post- 
collapse model. See Addendum C (section 5.8). 

The effects of anisotropy are studied in Fig. 5.7, where we can see that the external parts of the cluster 
are dominated by radial orbits. Inside the critical radius (indicated by a star symbol), one notices a 
slight tangential anisotropy, an effect of the depletion of loss-cone orbits. At large radii, the velocity 
distribution tends to isotropy as an effect of the outer bounding condition imposed at 10 4 pc. 
The effects of anisotropy in the stellar system can also be seen in Fig. 5.9, where we plot the components 
along the line of sight Olos (solid line) and on the sky ("proper motions"). The latter is decomposed 
into the radial direction (i.e. towards/away from the position of the cluster's centre) component, OpM.r 
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T= 9.63X10 1 T rh (0) = 1.00x10 ,0 yrs 
Radius in pc 
0.1 1 10 



~i — i — mm 




0.01 0.1 1 10 

Radius in Af-body units 

Fig. 5.7: Profile of the anisotropy parameter for the same case and time as Fig. 5.5. The decrease of at 
the border is an artefact of the inappropriate boundary condition. An outer boundary with radial anisotropy 
should be open, but here we enforce the adiabatic wall. If one "opens" the wall, but it would be at the 
expense of the stability of the program. All this does affect only a very small fraction of the total mass. 



(dashes) and the tangential component, OpM.t (dash-dot). Note that the radial anisotropy in the outskirts 
of the cluster reveals itself as a radial "proper motion" dispersion slightly larger than the other com- 
ponents. For an isotropic velocity dispersion, all three components would be equal. Despite loss-cone 
effects, there is no measurable anisotropy at the centre. 

The loss-cone induced anisotropy could be detected only if one could select the stars that are known to 
be spatially close to the centre (and not only in projection) as would be feasible these stars happened 
to be of a particular population. An interesting possibility that we'll soon investigate with multi-mass 
models is the concentration at the centre of more massive stars, i.e. mass segregation. 
Note that some anisotropy has been detected among the stars orbiting the central massive black hole of 
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T= 9.63X10 1 T rh (0) = 1.00x10 ,0 yrs 
Radius in pc 
0.1 1 10 




0.01 0.1 1 10 

Radius in Af-body units 

Fig. 5.8: Projected density for the same case and time as Fig. 5.5. In the interval between and r cr ; t we 
get a slope of —3/4, as expected. 



the Milky Way, Sgr A*, at distances closer than 1, i.e. 0.04 pc which is well inside the critical radius 
(> 1 pc) (Schodel et al., 2003). However, the detected anisotropy is in the radial direction rather than 
tangential. It is probably not connected to loss-cone effects but to particular history of these seemingly 
very young stars (Ghez et al., 2003) which remains a puzzle. 

In Fig. 5. 12 the diffusion model is layed out for the loss-cone as evaluated in previous sections. We eval- 
uate a model close to the post-collapse moment analysed in the other plots. We depict here the loss-cone 
filling factor K (upper panel), the loss-cone and diffusion angles 9 D and 0i c (middle panel) and the local 
contributions to the total loss-cone star accretion rate (lower panel). Our diffusion model reproduce well 
the picture of Frank and Rees (1976): the critical radius is defined by 9d = 0i c and coincides with the 
radius where the local contribution to the loss-cone accretion rate has its peak value. These two angles 
are connected to the time-scales f; n and t out , 9p « foutArelax and d\ c °< fin/frelax- The figures show that 
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Fig. 5.9: Projected velocity dispersions for the same case and time as Fig. 5.5. The line of sight component 
is represented with a solid line and the proper motion component with a dashed one for the radial and 
dashed-dot for the tangential contribution. 



the maximum contribution to the mass accretion rate stems at the radius where 9d — (he- Consistently, 
Frank and Rees (1976) estimated the total mass accretion rate as °< p(r cr j t ) rL t /f I ei a x('crit)- 
The dependence of loss-cone accretion rate on time during the late re-expansion phase can be estimated 
through simple scaling laws. One starts with the relation of Frank and Rees (1976) mentioned above, 

ss ^P( r crit)4t (5 35) 



?relax(?" ( 



ait) 
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Time in years 




1 10 10 2 

Time in T rh (0) 

Fig. 5.10: Evolution of the stellar density in the central region for our model with 10 5 stars and ^,h(0) = 
50M S (case a). The solid line depicts the density at the influence radius R^. The dashed line shows the 
average density within R^. 



and the definition of the critical radius, 

01c('"crit) = 0rj( r crit) 



M (5-36) 



r crit frelax (?"crit ) 

One substitutes the following relations into Eq. 5.36, 

r x oc J(l^ \ 



r 3/2 

Across ( r ) °^ 



?relax ( r ) ' 



1/2 ' (5.37) 

3 yy3/2 



n(r) r 3 / 2 n(r)' 
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Fig. 5.11: Same as Fig. 5.10 but for the three dimensional velocity dispersion. 



where we have made use of the fact that the potential is dominated by the BH in the region of interest 
(fait < Hi)- Finally one needs the dependence of the density of stars on time and radius, n(r, T). We have 
seen that, to a good approximation, the re-expansion of the cluster is homologous with Lagrange radii 
expanding like R °c r 2 / 3 . In the region between r cr i t and r^, the density profile resembles a power-law 
cusp; hence, a general self-similar evolution can be described by 

n(r,T)=no(T)(^-^) , (5.38) 

where ro is some Lagrange radius. Hence, from conservation of mass inside rn, 

n(r,T) oc T 2 ^ r -« . (5.39) 
Combining relations 5.36, 5.37 and 5.39, one finds 
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1 2(3-o) 



(5.40) 



and, inserting this into Eq. 5.35, 




(5.41) 



which, by integration, yields 




(5.42) 



— J\UL — 3) 



-s(q-3) 



(5.43) 



For a = 7/4, which is appropriate here (again because r cr ; t < r^), the exponent in the last relation turns 
out to be —95 /79 ~ — 1 .20, in remarkable agreement with figures 5.1b and 5.2b. 

5.5 Discussion 

We have presented in this chapter a method to follow the evolution of a spherical stellar cluster with 
a central accreting BH in a fully self-consistent manner concerning the spatial resolution. As regards 
the velocity space, we use a simplified model based on ideas of Frank and Rees (1976) in order to 
describe the behaviour of the distribution function inside and outside the loss-cone by a simple diffusion 
equation. This numerical method is an extension of the "gaseous model" which has been successfully 
applied to a variety of aspects of the evolution of globular clusters without central BH (Spurzem and 
Takahashi, 1995; Spurzem and Aarseth, 1996; Giersz and Spurzem, 2000, 2003; Deiters and Spurzem, 
2001, amongst others). With this new version, the simulation of galactic nuclei is also feasible. 
In addition to an explanation of the physical and numerical principles underlying our approach, we have 
concentrated on a few simple test computations, aimed at checking the proper behaviour of the code. 
We considered a system where all stars are and remain single, have the same mass, stellar evolution 
and collisions are neglected and a seed central BH is allowed to grow by accreting stellar matter through 
tidal disruptions. The present version of the code already allows for a (discretised) stellar mass spectrum 
and stellar evolution and we are in the process of including stellar collisions because they are thought 
to dominate over tidal disruption in most galactic nuclei, as far as accretion on to the BH is concerned 
(David et al., 1987a,b; Murphy et al., 1991; Freitag and Benz, 2002). In a subsequent chapter, we shall 
increase complexity and realism one step further and consider systems with a mass spectrum. Using 
both this gaseous code and the Monte Carlo algorithm (Freitag and Benz, 2001, 2002), we will investi- 
gate the role of mass segregation around a massive black hole (Amaro-Seoane, Freitag & Spurzem, in 
preparation), a mechanism which may have important observational consequences as it probably affects 
the structure of the central cluster of the Milky Way (Morris, 1993; Miralda-Escude and Gould, 2000; 
Freitag, 2003b, a; Pfahl and Loeb, 2003) and impacts rates of tidal disruptions and capture of compact 
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T = 9.95x 1 0' 7 rh (0) = 1 .04x 1 10 yrs 

log, R [pc] 
-10 12 




log 10 ff [N-body] 

Fig. 5.12: In this triple diagram we evince the dependence on radius of the mass accretion rate, the loss-cone 
and diffusion angle which are related to the filling and depletion time-scales of the loss-cone (see text) and 
the filling factor K. The critical radius is defined by the condition 0o = 8\ c (broken line). 

stars by emission of gravitational waves in dense galactic nuclei (Magorrian and Tremaine, 1999; Syer 
and Ulmer, 1999; Sigurdsson, 2003, and references therein). 

Unfortunately, the literature has relatively little to offer to check our models. The most robust predictions 
are probably the analytical and semi-analytical analysis for the regime where the gravity of the BH 
dominates and the Fokker-Planck treatment of relaxation holds (Bahcall and Wolf, 1976; Shapiro and 
Lightman, 1976; Lightman and Shapiro, 1977; Cohn and Kulsrud, 1978). The most important feature 
of these solutions is that, provided the system is well relaxed and one stands beyond the critical radius 
(inside of which loss-cone effects complicate the picture), a cuspy density distribution is established, 
p«r a with a = 7/4. Our code nicely agrees with this prediction. 

Concerning the evolution of the system, we first note that, initially, the cluster follows the usual and 
well understood route to core-collapse. That the gaseous model can successfully simulate this phase has 
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been clearly established in previous works (Giersz and Spurzem, 1994; Spurzem and Aarseth, 1996). 
When the core has become dense enough, the BH starts growing quite suddenly. As it accretes stars that 
are deeply bound, i.e. with very negative energies, the BH creates a outward flux of energy and allows 
the cluster to re-expand. As long as the source of energy is centrally concentrated and that the mass 
of the BH remains relatively slow, one expect the re-expansion to become self-similar, a regime during 
which the size of the cluster increases like R °c T 2 / 3 (Henon, 1965; Shapiro, 1977; McMillan et al., 
1981; Goodman, 1984, among others). This is again well reproduced by the gaseous model. Solving 
the Fokker-Planck equation with a Monte Carlo method, Marchant and Shapiro (1980) and Duncan and 
Shapiro (1982) have realised a series of simulations of single-mass globular clusters with a central BH. 
Because their resolution was quite low and because they used "initial" conditions difficult to implement 
(in most of their runs the central BH is not present initially but introduced at some instant during deep 
collapse), we do not attempt a quantitative comparison with their results. An added difficulty is that we 
do not include tidal truncation of the cluster. However, an important finding of Marchant and Shapiro 
(1980) is reproduced by our computations, namely that the initial mass of the seed black hole has little 
effect on the post-collapse evolution, provided it represents only a small fraction of cluster mass. In 
particular, the BH mass at late times converges to the same value which only depends on the size and 
mass of the cluster. We note that such convergence was also obtained with the Monte Carlo algorithm 
and that a comparison between results obtained with that code and an early version of the program 
described here was presented in Freitag and Benz (2002). More comparisons between the two methods 
are planned (Amaro-Seoane, Freitag & Spurzem, in preparation). 

Here one has to mention that the energy input of the BH star accretion causes a "temperature" increase in 
the central region which is followed by a thermal expansion. Therefore, the system is a normal thermal 
system with positive specific heat in contrast to the cores of self-gravitating systems, where the energy 
input due to binary hardening causes a core expansion and decrease of temperature (Bettwieser and Sug- 
imoto, 1984). Afterwards, an inversion in the radial temperature profile follows and the expansion is a 
reverse gravothermal instability (Spurzem, 1991). Since our system is dominated by external gravitation 
in the centre we cannot expect such a behaviour. Furthermore, as the central BH grows irreversibly, it 
continues accreting stars in spite of the re-expansion and continuous decrease of the density of stars. 
Hence, a second core collapse is impossible and oscillations of the central cluster density do not occur. 
Among the aspects of our results that require further investigation, we mention the shape of the density 
profile inside the critical radius. Although the well known p r~ 7 / 4 Bahcall-Wolf solution only strictly 
applies for r cr i t < r < r\, this cusp extends inward nearly down to the tidal disruption radius in the 
stationary models of Marchant and Shapiro (1979), which also include loss-cone physics. This result 
is in disagreement with the analysis of Dokuchaev and Ozernoi (1977) (see also Ozernoi and Reinhardt 
1978) which predicted p ^ r - l l 2 in this region. As shown on Fig. 5.5, we obtain an even stronger 
flattening of the density law inside r cr jt- At the present time it is unknown to us which solution, if any, is 
the correct one. A possibility to be considered is that this a consequence of the truncation of the moment 
equations to the second order. In other words, in regions where the loss-cone is significantly depleted, 
representing the velocity distribution by a simple dispersion ellipsoid and using the velocity dispersion 
to determine an "effective" loss-cone aperture (Eq. 4.85) is clearly quite a strong approximation. This 



120 



5.6 Addendum A: MBH wandering in a cuspy cluster. 



may impact the density distribution as the system adjusts its central structure to produce the heating rate 
required by the overall expansion. 

Fortunately, numerical inaccuracies at very small radii are unlikely to affect the overall structure and 
evolution of the cluster because the loss-cone accretion physics are essentially determined by the condi- 
tions at the critical radius, not in the immediate vicinity of the BH. 



5.6 Addendum A: MBH wandering in a cuspy cluster. 

Here we present a simple estimate of the wandering radius R wan of a MBH embeded in a stellar cluster 
whose density posseses a power-law cusp in the inner regions. We assume that, were it not for the effect 
of the MBH itself, the stellar cluster would be described by an eta-model (Dehnen, 1993; Tremaine 
et al., 1994) with enclosed stellar mass 

JZ*{r)=JtJ-^-X , (5.44) 
\l + r/aj 

where ^# c ] is the total mass in stars and a the break radius. For r <C a, the density is p « r~ a with 
a = 3 — 77. Inside /? W an, the MBH strongly perturbates stellar orbits and we suppose the density is 
rendered more or less constant. Hence, the potential felt by the MBH is approximately harmonic, 

*(r) = *o + J©V , with co 2 = G ^f wan) . (5.45) 

^ "wan 

For a harmonic oscillator, the RMS amplitude of the oscillations in velocity and space are linked to each 
other, 

V 2 - m 2 R 2 ~mfi> 2 - G-^* (ft wan) 

^RMS - 03 "RMS ~ G^wan ~ n • (5.40) 

"wan 

Dorband et al. (2003) have verified with N-body simulations that equipartition of kinetic energy between 
the MBH and the stars is established, at least in the case of 77 = 1.5. Namely, 

Jt.V^-m^o 2 , (5.47) 

where a is the stellar velocity dispersion at r — a. For 77 = 1.5 (Tremaine et al., 1994), 

a 2^ 01 G^el 
a 

Finally, assuming /? W an <C a and, hence, .-#*(7?wan) ~ ^# c i ftwan/a) 1 ' and combining equations 5.46, 
5.47 and 5.48, we obtain 

fl wan w 0.01a f J for 77 = 1.5 (5.49) 
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and 

U) (5.50) 

for general eta-models. 



5.7 Addendum B: Velocity dispersion in the central regions 

In the region dominated by the central BH, one may expect a "Keplerian" profile for the velocity dis- 
persion, a <x r~ l l 2 . However, in Sec. 5.4, we have seen that in our standard model (10 5 stars, case a this 
relation does not really apply where expected, i.e. between the "1-star" and the influence radius. Here 
we show that the spherical Jeans equation for a system in dynamical equilibrium is nevertheless obeyed. 
In spherical symmetry the assumption of dynamical equilibrium (or stationarity) amounts to u = (v r ) = 
0. We note that the gaseous model can cope with u ^ 0. On the other hand, for a system whose evolution 
is driven by relaxation, one expects a r ,t S> u. The Jeans equation then reads (Binney and Tremaine 1987, 
Eq.4-55) 

? /dlnn din a? „„\ 
GM r ^-ra 2 -— + — -^+2/3 . (5.51) 
\dlnr dlnr / 

M r is the mass enclosed by the radius r, n the number density of stars and 2j3 = 2 — a 2 jo 2 is the 
anisotropy parameter; other quantities have been defined previously. One sees easily that if M r = Ji* 
and the first and third term in the brackets are both constant, a r~ l l 2 at small r. But, as figures 5.13 and 
5.14 demonstrate, none of these assumptions exactly applies in the range of radius under consideration. 
Consequently, we do not get a clean Keplerian velocity profile although (or because) Eq. 5.5 1 is satisfied. 
Finally, we mention that our models with 10 6 and 10 7 stars (and same initial size) exhibit a Keplerian 
velocity cusp outside the 1-star radius, during their post-collapse evolution. This is partly due to the 
relatively more massive black hole (larger influence radius) and partly to the much smaller 1-star radius. 



5.8 Addendum C: Projected velocity dispersions 

For the Figures 5.9 and 5.8 we have integrated the density along the z-axis for the projection, 

£(r)= / p(Vr 2 + Z 2 )dz = 
Jz=0 

2 l p{R) vw^ dR (5 ' 52) 

If one observes the cluster along the z-axis, the contributions to the projected velocity dispersions are 
as indicated in Fig. 5.15. We have that R = \/r 2 +z 2 and Gg = o§ = a t , where the subscript t stands for 
tangential. We can reckon that 
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Fig. 5.13: Check of stationary Jeans equation for our "standard" model (10 5 stars, case a) at T = WGyrs 
(same model and time as Fig. 5.6). Panel (a) depicts the logarithmic derivatives of the stellar density n and 
radial component of the velocity dispersion well as the anisotropy parameter, 2/3 = 2 — a} /a}. Horizontal 
lines corresponding to values of 0, —1, —1.75 and —2.23 are present to guide the eye. In particular, the 
"Keplerian" velocity profile is din cr^/dln r = —I, The round dot indicates the influence radius, the star the 
critical radius and the square the "1-star" radius. On panel (b), we plot the three terms of the right side of 
the stationary Jeans equation and check that their sum is (nearly) equal to the left side term, i.e. GM r . 
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Fig. 5.14: Blow-up of panel (b) of Fig. 5.13 for the region between the 1-star and the influence radius. 



°z 2 = °r cos2 + °f s i n2 $ 

a r 2 = CRsin 2 e + a t 2 cos 2 e, (5.53) 

where we have defined a 2 = (J t 2 /2, to be consistent with the notation used until now. 
In Fig. 5.15 we can see that (7 Z contributes to <7los> o r to Op m , r and o 2 to Cp m t. 
Thus, we obtain the projected velocity dispersions, 



OlosW = J z=Q (ol(R) cos 2 9 + a? (R) sin 2 9) p(R)dz = 



2 



K/?)i=o (|i( C7 rW-^W) + ^ 2 W)pW* 
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2 fRmdx 

a^W = ^ =o (OrW sin 2 + & t 2 (/?) cos 2 d)p (R)dz = 

" -y z _ (^(& t 2 w-^w)+^(fi))p(«)* 



a ^ {r)= W)Lo dt (i?)p(fi) *' (5 - 54) 

since sin 2 = r 2 //? 2 = 1 - cos 2 9 and cos 2 = z 2 /R 2 . 
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Chapter 6 



Multi-components clusters with/out a 
central BH: Mass segregation 

6.1 An academic exercise: Mass segregation in two mass-component 
clusters 

SINCE we have arduously examined the idealised situation for a stellar cluster in which all stars 
possess the same mass in order to subject to scrutiny our model, we have morally the right to 
extend the analysis a further step. Here we tackle with more realistic configurations in which 
the stellar system is splitted into various components 1 . The second integer number immediately after 
one is two and so will we first extend, cautious and wary as we are, our models to two-components 
star clusters. The initial motivation for the study of such physical systems was the better resemblance 
to real clusters. The processes that one-component clusters bring about is nowadays well understood 
and has been plenteously studied by different authors to check for the goodness of their approaches. 
New features of these systems' behaviour arise when we consider a stellar system in which masses are 
divided into two groups. Depending on how the system taken into consideration is configurated will we 
exclude dynamical equilibrium (here it is meant that the system is not stable on dynamical time-scales) 
or equipartition of different components kinetic energies is not allowed (thermal equilibrium). 
Spitzer (1969) gave the analysis an initial shove with his study. For some clusters it seemed impos- 
sible to find a configuration in which they enjoy dynamical and thermal equilibrium altogether. The 

'This chapter is part of present work and will be published like Freitag, Amaro-Seoane and Spurzem (2004). Some part of it 
has been employed for Khalisi, Amaro-Seoane and Spurzem (2004) -submitted to MNRAS- and an initial concept was employed 
for Amaro-Seoane and Spurzem (2003) 
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heavy component sink into the centre because they cede kinetic energy to the light one when reaching 
equipartition. The process will carry on until equipartition is fully gained on and the heavy component 
self-gravity moves out in the core. In the most of the cases, equipartition happens to be impossible, be- 
cause the subsystem of massive stars will undergo core collapse before equipartiton is reached. Anon, a 
gravothermal collapse will happen upon this component and, as a result, a small dense core of heavies is 
formed (Spitzer, 1969; Lightman and Fall, 1978). This gravothermal contraction is a product of negative 
heat capacity, a typical property of gravitationally bound systems (Elson et al., 1987). 
Disparate authors have addressed the problem of thermal and dynamical equilibrium in such systems, 
from recent direct jV-body simulations (Portegies Zwart and McMillan, 2000) and Monte Carlo simu- 
lations (Watters et al., 2000) to direct integration of the Fokker-Planck equation (Inagaki and Wiyanto, 
1984; Kim et al., 1998) including Monte Carlo approaches to the numerical integration of this equation 
(Spitzer and Hart, 1971). For a general and complete overview of the historical evolution of two-stars 
stellar components, see Watters et al. (2000) and references therein. 

If we do not have any energy source in the cluster and stars do not collide (physically), the contraction 
carries on self-similarly indefinitely; in such a case, one says that the system undergoes core-collapse. 
This phenomenon has been observed in a big number of works using different methods (Henon, 1973; 
Henon, 1975; Spitzer and Shull, 1975; Cohn, 1980; Marchant and Shapiro, 1980; Stodolkiewicz, 1982; 
Takahashi, 1993; Giersz and Heggie, 1994; Takahashi, 1995; Spurzem and Aarseth, 1996; Makino, 
1996; Quinlan, 1996; Drukier et al., 1999; Joshi et al., 2000, etc) Core collapse is not just a characteristic 
of multi-mass systems, but has been also observed in single mass analysis. 

Spitzer (1969) gave the analytical criterion to determine whether a two-component system has achieved 
energy equipartition. According to his analysis, energy equipartition between the light and heavy com- 
ponent exists if 



Where and are the total numbers of light and heavy components, respectively. More recent 
numerical calculations (Watters et al., 2000) have settled this criterion to A, 



When we modify the ratio ^# m ax/^, the time required to reach core-collapse is different. In a cluster 
with, for instance, a broad Salpeter IMF between [0.2M©, 120 M ] core-collapse takes place after < 
0.1 frh(0), whereas for a single-mass Plummer model it occurs after > 10f r h(0) (this example was taken 
from the Monte Carlo-based calculations of Giirkan et al. 2004) 

There is an ample evidence for mass-segregation in observed clusters. McCaughrean and Stauffer (1994) 
Hillenbrand and Hartmann (1998) provided a new deep infrared observations of the Trapezium cluster in 
Orion that clearly show the mass segregation in the system, with the highest mass stars segregated into 
the centre of the cluster. To test whether there is evidence for more general mass segregation, they show 
in Fig. (6.1) cumulative distributions with radius of stars contained within different mass intervals. They 
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Fig. 6.1: In this plot by Hillenbrand and Hartmann (1998) we have a clear-cut evidence for mass-segregation 
of stars more massive than 5 Mq (long-dashed lines) toward the cluster centre and some evidence for general 
mass segregation persisting down to 1-2 M Q in the Orion Nebula cluster. The cumulative radial distributions 
of source counts over different mass intervals are shown. To clarify the sensitivity of the cumulative plots 
to the outer radius they have shown here four panels with four different limiting radii 



include in the plot four different panels in order to manifest the sensitivity to the limiting radius. They 
find that, inside 1 .0 pc, general mass segregation appears to be established in the cluster, with stars of 
masses less than 0.3, 0.3-1.0, 1.0-5.0, and greater than 5 M Q progressively more centrally concentrated 
with increasing mass. 

At this point, the question looms up whether for very young clusters mass segregation is due to relax- 
ation, like in our models, or rather it reflects the fact that massive stars are formed preferentially towards 
the centre of the cluster, as some models predict. 

Raboud and Mermilliod (1998) study the radial structure of Praesepe and of the very young open cluster 
NGC 623 1 . There they find evidence for mass segregation among the cluster members and between 
binaries and single stars. They put it down to the greater average mass of the multiple systems. In 
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Fig. 6.2: Mass segregation in NGC 623 for two mass interval sets. The two left panels include all sample 
stars, whereas the right ones do not include the 9 bright stars of the cluster corona. For the two top 
figures M < 5M S (filled squares), M € [5, 1O[M (open squares), M e [10, 2O[M (crosses) and M > 20M S 
(triangles). For the two bottom figures, M < 2.5M Q (filled squares), Me [2.5, 6.3[M© (open squares), 
Me [6.3, 15.8[M (crosses) and M> 15.8M Q (triangles) 



Fig. (6.2) of Raboud and Mermilliod (1998) again we have clear evidence for mass segregation in NGC 
623 1 . In the two first panels the mass intervals are set in a different way to those in the bottom. 
The two left-hand panels include the 9 bright stars of the cluster Corona, while the right hand two 
diagrams do not. The manifestation of mass segregation for massive stars (triangles) is clearly nailed 
down, whereas for stars with masses e [5, 20] M are spatially well mixed (open squares and crosses); 
i.e., mass segregation is not yet established over a rather large mass interval. This population is more 
concentrated than the lower-mass population (here shown with filled squares). They make out from this 
Fig. (6.2) that only a dozen, bright, massive, mainly binary stars are well concentrated toward the cluster 
centre. 
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It seems therefore interesting to set out for multi-mass models with the two-component ones as a start- 
ing point to take care of, since it is well-studied and we have robust observational proofs of this phe- 
nomenon. On the other hand, observations do not tell us whether mass segregation is due to relaxation. 
For this aim, and as the title of this section indicates, we have performed here a somehow homiletic 
but in any case interesting exercise in order to make an acid test of the gaseous model. In less than 
three days we performed a whole set of 10 4 simulations for two-component models. We define two 
parameters now that describe the physics of the system, 

q := -M^j 'M ', H :— m^/m\ (6.3) 

In this definition, ^ is the total mass of the system, Jt^ the total mass in heavies and mh,i the mass of 
one heavy (light) star. 

If £ = 1 — q, we have let £ vary from 10~ 4 to 9.99 • 10 -1 . For each £ value, we let jj, vary between 1.03 
and 10 3 . 

The values for q are regularly distributed in log(£). For £ w 1 we have added a series of values in 
log (£ — 1). The mean particle mass is 1M and the total mass 1O 6 M , but this has not importance for 
our study, because the physics of the system is ruled by relaxation and therefore the only relevant thing 
is the relaxation time. If we are using Fokker-Planck or half -relaxation units, we can always extend the 
physics to any other system containing more particles (if only relaxation is on play). As a matter of fact, 
one relates the N-body units to Fokker-Planck units as follows: 

^FP = %B-<A/ln(r-^C) (6-4) 
Therefore, the relation between FP- and jV-body- time units is 

^Fpltime - %B|time-<A/ln(y-^) (6.5) 

In order to be able to compare our results with N-body, this relation is very useful. The mean mass is 
therefore just a normalisation. What really determines the dynamics of the system are the mass ratios, q 
and jj,. 

In Fig. (6.3) we show the whole (q, /x)-parameter space in a plot where the time at which the core- 
collapse begins is also included. The green zone corresponds to the quasi single-mass case. In the red 
zone we have the largest difference between masses and blue is an intermediate case. 
In Fig. (6.4) we show collapse times for clusters models with two mass components normalised to the 
single-mass core-collapse time r cc (s.s.) for different values of ji. The initial clusters are Plummer 
spheres without segregation. The collapse times are displayed as a function of the mass fraction of the 
heavy component in the cluster. When compared to single-mass component systems, we see that the 
core-collapse time is accelerated notably for a wide range of the heavy component Ji^ (Mi). Even a 
small number of heavies accelerate the core-collapse time. 

It is really interesting to compare the capacity of our approach by comparing the results of this set of 
simulations to the N-body calculations of star clusters with two-mass components performed by Khalisi 
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T cc (FP time) 

■ 4.70e-04 - 1,33e-03 

■ l.34e-03- 2.98e-03 

■ 2.98e-03 - 5.84e-03 

■ 5.85e-03 - 1 .05e-02 

■ 1.Q5e-02- 1.76e-02 

■ l.77e-02 - 2.84e-02 

■ 2.84e-02 - 4.48e-02 

■ 4.48e-02 - 6.80e-02 

■ 6.81e-02 - 1.01e-01 

■ l.01e-01 - 1.45e-01 

■ 1.45e-01 - 2.02e-Q1 

■ 2.02e-01 - 2.73e-01 

■ 2.74e-01 - 3.57e-01 

■ 3.58e-01 - 4.53e-01 

■ 4.54e-01 - 5.59e-01 

■ 5.59e-01 - 6.66e-01 

■ 6.67e-01 - 7.78e-01 

■ 7.78e-01 -S.88e-01 

■ 8.88e-01 -9.91e-01 

■ 9.91e-01 - 1.09e+00 

■ l.09e+00- 1.17e+00 

■ l,17e + 00- 1.24e+00 

■ L24e+00- 1.30e+00 

■ l.30e+00- 1.35e+00 

■ L35e+00- 1.39e+00 

■ 1. 3964-00- 1.41e+00 
BL41e+00- 1.43e+00 

■ l.43e+00 - 1.44e+00 

■ l.44e + 00- 1.45e+00 

■ 1 .45e+00 - 1 .46e+00 

12 3 

log(^) = log(m h /m 1 ) 

Fig. 6.3: Parameter space for the set of 10 4 simulations. Here t en & stands for the core collapse time and is 
expressed in FP units (see text); time at which the simulation ended, q and fl are plotted logarithmically. 



(2002) during the completition of his doctoral thesis. For this aim, we plot the evolution of the average 
mass in Lagrangian shells of the cluster from the averaged mass in Lagrangian spheres containing the 
following mass percentages [0 — 1] , [2 — 5] , [10 — 20] , [40 — 50] , [75 — 95] %, among others, to be able 
to compare with the results of Khalisi. These are the comprised volume between two Lagrangian radii, 
which contain a fixed mass fraction of the bound stars in the system. 

We have calculated the average mass as follows: If is the total mass for the component i comprised 
at the radius r and is the average mass for this component within that radius, we can find out what is 
the value of m* 1+1 ) (the average mass between m* and m£ 1+ 1 ^ ) knowing Af, , and m£ +1 \ 
This is schematically shown in Fig. (6.5). Indeed, 

M r (i+1) = NP ■ + • +1) = Jvf +1) • (6.6) 

Since 
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Fig. 6.4: Core-collapse time for different values of q and fl 



where 



we have that, from Eq. (6.6), 



A^ (l+1 W r (i) +Af (6.7) 



NP = *^ (6.8) 



M r ^-M« (69) 



We show in Figs. (6.6) and (6.7) we show the curves corresponding to the values shown in table 6.1. 
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Fig. 6.5: Average mass in Lagrangian shells from averaged mass in Lagrangian spheres 



jU in Khalisi (2002) 


jU in this work 


1.25 


1.27 


1.5 


1.56 


2 


2.06 


3 


2.92 


5 


5.09 


10 


10.2 



Table 6.1: Different fl values used in the iV-body calculations and in our models for Fig. (6.6) 

We have followed in the curves the evolution of the system until a deep collapse of the system. They 
show the evolution until the most massive component dominates the centre. 

In order to compare our plots with those of Khalisi (2002), one should look in his diagrams in the region 
during core contraction. At this point, we can observe in Fig. (6.6) a self-similarity after core-collapse 
(Giersz and Heggie, 1996). Binaries are responsible for interrupting core-collapse and driving core re- 
expansion in the A^-body simulations. The flattening in the A^-body plots at the moment of core-collapse 
is due to the binary energy generation. This means that we can only compare the steep rise, but not the 
saturation. 

For instance, in the second plot of the A^-body set (second column on the top), we have to look at the 
point in which the average mass of the N-body system is about 1 .20 in the — 1 % shell. This establishes 
the limit until which we can really compare the behaviour as given by both methods. Our simulations 
yield a very similar evolution until that point. The gaseous model behaves (it clearly shows the tendency) 
like Af-body. 

From Eq. (6.5), we get that the conversion factor is the same; namely, for y = 0.11, \\\{y ■ / -^K = 
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B: Mean star moss in shells 
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Fig. 6.6: Average Lagrangian radii shells for the A'-body models of Khalisi (2002) (see text for further 
explanation) 
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Fig. 6.7: Average Lagrangian radii shells for our models, equivalents to those of Fig. (6.6) 
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0.0022. On the other hand, the value of y is not so well defined and depends on the cluster's mass 
spectrum (Henon, 1975). This means that potentially it is not the same for the different models. For a 
broader mass spectrum, y is about 0.0 1 and, unfortunately, in the case of having a small particle number, 
it will definitively make an important difference despite the "smoothing" effect of the logarithm, viz 
ln(y- jY-t)l jY* = 0.0013. Thus, in order to be able to compare the different models, one should consider 
y as a free parameter ranging between 0.01 to 0.2 and look for the best fit for the most of the cases. On 
the other hand, we must bear in mind that the jV-body simulations of Khalisi (2002) do not go in deep 
core collapse and, so, the moment at which the core radius reaches a minimum is not the same as for 
our model. To sum up, we cannot really exactly say until which point we can compare, because the core 
collapse time will be different. 

6.2 Clusters with a broader (> 2) mass spectrum without a BH 

So as to be able to interpret observations of young stellar clusters extending to a larger number of mass 
components, it is of paramount relevance to understand the physics behind clusters without a central BH 
first. It has been shown that for a cluster with a realistic IMF, equipartition cannot be reached, for the 
most massive stars build a subsystem in the cluster's centre as the process of segregation goes on thanks 
to the kinetic energy transfer to the light mass components until the cluster undergoes core collapse 
(Spitzer, 1969; Inagaki and Wiyanto, 1984; Inagaki and Saslaw, 1985). 

Whereas the case in which the BH ensconces itself at the centre of the host cluster is more attractive 
from the dynamical point of view, one should study, in a first step, more simple models. 
In this section we want, thus, to go a step further and evaluate stellar clusters with a broad mass function 
(MF hereafter). We study those clusters for which the relaxation time is short enough, because this will 
lead the most massive stars to the centre of the system due to mass segregation before they have time to 
leave the main sequence (MS). In this scenario, we can consider that stellar evolution plays no role; stars 
did not have time to start evolving. The configuration is similar to that of Giirkan et al. (2004), but they 
employ a rather different approach based on a Monte Carlo code, using the ideas of Henon (1973) that 
allows one to study various aspects of the stellar dynamics of a dense stellar cluster with our without a 
central MBH. Our scheme, although being less realistic than MC codes (and A^-body ones) and unable, 
in its present version, to account for collision has the advantage, as already mentioned in fore-going 
sections, of being much faster to run, and providing data that has no numerical noise. It captures the 
essential features of the physical systems considered in our analysis. 

6.2.1 Mass segregation in realistic clusters 

One of the first questions we should dope out is up to how many components one should take into 
consideration when performing our calculations. Since the computational time becomes larger and 
larger when adding more and more components to the system, we should first find out what is a realistic 
number of components in our case. For this end we performed different computations with different 
number of stellar components. 
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For the simulations shown here, the initial cluster models are Plummer models with a Salpeter MF 



between 0.2 and 120 M . In this equation a takes a = 2.35 (Salpeter). There is no initial mass segre- 
gation. 

The discretisation of the mass components has been done as follows: 



In this equation 8 is the discretisation exponent. If 8 > 1 we have more bins at low mass; for 8 < 
1, we have more bins at high mass. I.e., 8 allows one to put more discretised mass components at 
low masses (8 > 1) or at high masses (8 < 1), 5=1 gives the logarithmical equal spacing. M max m i n 
are, respectively, the maximum and minimum individual stellar masses for the components. For all 
simulations the number of mass bins has been typically set to 15. We have chosen a Plummer model 
by default and the models have 10 6 stars. The model radius by default is Rpi = 1 pc. The default initial 
mass function is a Salpeter. 

In Fig. (6.8) we show the Lagrangian radii for ten different models and look for main dynamical char- 
acteristics of the system's behaviour, the core collapse time, the Lagrangian radii containing 90, 70, 50, 
20, 10, 3, 1, 0.3, 0.1, 3 • 10~ 2 , 10~ 2 , 3 • 10~ 3 and 10~ 3 % of the stellar mass. 

In this plot, A^ C omp stands for the mass spectrum different components number. For A^ om p = 6 we have 
performed three simulations varying the 8 parameter between 1 .0 (equal logarithmic spacing of compo- 
nents), 0.75 (more massive components) and 0.5 (even more). For N com p = 12 we have performed one 
only simulation (with 8 = 1, by default); for the A^omp = 20 case we have repeated the same procedure 
as with six components, the last but one that we have chosen is N comp = 20 and, in this case, we stud- 
ied two grid resolutions, A^h = 200 (the default value) and 400 grid points, in order to check whether 
this could influence the results (see chapter 2). To finish with, a big simulation with Acomp = 50 was 
performed and included in the analysis. Whereas we can see an important difference between models 
of 6 and 12 components, we see that the global behaviour from 12 components onwards is very simi- 
lar. Therefore, unless indicated, we choose 15 components in our study in this section, since a higher 
number would not contribute anything essential. 

To see this in more detail, in Fig. (6.9) we show the Lagrangian radii for each stellar mass m, and the 
corresponding mass fraction f m for the 25 and 15 components simulations. Again, we cannot see any 
substantial difference between the 25 and 15 cases. 

Taking the last arguments into account, we have done an analysis of mass segregation in multi-mass 
models with more than two stellar components without BH. In Figs. (6.10) and (6.1 1), we show in the 
evolution of a stellar cluster of 15 components (in colours); m is the mass (in M ) of the stars in each 
component and f m the corresponding fraction of the total mass. In the upper box we have the density 
profile, where the solid black line represents the total density; below, we have the average total mass for 
the system. We show different moments of the system. At T = we have the initial model, which duly 




M, 



<5 



(6.11) 
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Fig. 6.8: Lagrangian radii and average stellar mass for 10 models with different mass spectrum (see text) 
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Fig. 6.9: Lagrangian radii for each stellar mass m, in the cases of 25 and 15 mass components 
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T = 0.00e + 00T rh (0) 
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Fig. 6.10: Initial density profile for a stellar cluster with 15 components (upper panel) in ,/V-body units and 
average total stellar mass in Mq (lower panel) 



shows no mass segregation. As time elapses, at T = 5.30 ■ 10~ 2 7i-h(0), we observe how mass segregation 
has fragmented the initial configuration; the heavy components have sunken into the central regions of 
the stellar cluster and, thus, risen the mean average mass. The outer parts of the system start losing 
their heavy stars quickly and, consequently, their density profile retrogresses. This becomes more acute 
for later times at T = 6.75 • 10 _2 7;-h(0), as the right plots of Fig. (6.1 1) show. In these plots and, more 
markedly in the right panel of density profile, we can observe a depletion at intermediate radius. 

6.2.2 Core-collapse evolution 

Gtirkan et al. (2004) shown that for a broad MF -which can be a Salpeter or a Kroupa-, mass segre- 
gation produces a core-collapse of the system that happens very fast. For clusters of moderate initial 
concentration, they find that this is in about 10 % of the Ti-h(O), the initial half-mass relaxation time (i.e. 
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T = 6.75e-02T rh (0) 




Fig. 6.11: Same situation as in Fig. (6.10) but at later times. See text for explanations. 
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Fig. 6.12: Illustration of core-collapse in multi-mass systems treated with a Monte Carlo approach (courtesy 
of M. Freitag) 

the half-relaxation time that the cluster had when time started, at t = 0). A good and clear illustration of 
this is Fig. (6.12) and Fig. (6.13). In the former one, on the left panel we have the initial configuration 
of the system. On the right one, we have the cluster at the moment of core-collapse. In the figure, all 
stars within a slice containing the centre has been depicted. On the other hand, this does not represent 
a real physical system, because all stars' radii have been magnified (see bottom of each panel). The 
dashed circles represent spheres containing 1, 3 and 10% of the total cluster mass (from the centre). 
We can clearly see how the massive, large stars are segregated towards the centre. In Fig. (6.13), we 
show the core-collapse evolution of a multi-mass stellar cluster simulated with the gaseous model. As 
usual, m is the mass (in Mq) of the stars in each component and /„, the corresponding fraction of the 
total mass. On the left panel we display the time evolution of the central density for a model in which 
we have employed 15 individual mass components. The total density is given by the dotted line. On the 
right panel we have the evolution of the central velocity dispersions. The dotted black line shows the 
mass-averaged value 



jV-body units are used for the y-axes. 

One notes that, during core collapse, the central regions of the cluster become completely dominated by 
the most massive stars. But, contrary to the case of single-mass clusters, the central velocity dispersion 
decreases 2 (see Fig. 6.13). 

2 Priv. comm. M. Freitag 




m 



(6.12) 
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Fig. 6.13: Evolution of the central density and 3£>-velocity dispersion in a model with 15 components (see 
text for further explanations) 



6.3 Clusters with a broader mass spectrum with a BH 

The logical next step to the systems studied in previous sections is that of a multi-mass component 
cluster harbouring a central seed BH that grows due to stellar accretion. 

In this final section we extend our analysis to systems for which we use an evolved mass function of an 
age of about 10 Gyr. We consider a mass spectrum with stellar remnants. We employ a Kroupa IMF 
(Kroupa et al., 1993; Kroupa, 2001) with Mzams 3 from 0.1 to 120 M with the turn-off mass of 1M . 
We have chosen the following values for the exponent according to the mass interval, 

f 1.3, 0.008 < m*/M < 0.5 
a = I 2.2, 0.5 < m*/M < 1 (6.13) 
[ 2.7, 1 <m*/M < 120. 

And with the following component's distribution, 

• Main sequence stars of 0.1 — 1 M (~ 7 components) 

• White dwarves of ~ 0.6 M (1 component) 

• Neutron stars of ~ 1.4M (1 component) 

• Stellar black holes of ~ 1OM (1 component) 

The defined IMF evolves and gives us at later times an evolved populations with compact remnants. 
This means that main sequence stars can transform into white dwarves, neutron stars or stellar black 
holes according to their masses. If m^s is the mass of a MS star, we have defined the following mass 
ranges for the evolution into compact remnants: 

3 The zero age main sequence (ZAMS) corresponds to the position of stars in the Hertzsprung-Russell diagram where stars 
begin hydrogen fusion. 
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White dwarves in the range of 1 < m MS /M < 8 
Neutron stars for masses 8 < m M s /M Q < 30 
Stellar black holes for bigger masses, > 3OM 

As we already mentioned, we put at the centre a seed BH whose initial mass is of 50 M©. The initial 
model for the cluster is a Plummer sphere with a Plummer radius R P \ = 1 pc. The total number of stars 
in the system is = 10 6 . 

The presence of a small fraction of stellar remnants may greatly affect the evolution of the cluster and 
growth of the BH because they segregate to the centre from which they expel MS stars but, being 
compact cannot be tidally disrupted. This kind of evolution is shown in Figs. (6.14) and (6.15). 
Fig. (6. 14) shows us the time evolution of different Lagrange radii with 0.1, 10, 50, 80 % of the mass of 
each component. Here the core collapse happens at about T = 0.18r r h(0). The henceforth re-opening 
out is due BH accretion. 

In Fig. (6.15), we plot the density profiles of the system before and after the post-collapse phase. We 
can see that the slope of p « R^ 1 / 4 on account of the cusp of stellar BHs that has formed around the 
central BH. We can see how the different components redistribute in the process. We can see how the 
BH dominates the dynamics at the centre. 

We can study how the system evolves from the point of view of the distribution of kinetic energies 
between the different components of the clusters during the process of mass segregation. 
In Fig. (6.16) we show the evolution of the "temperature" of the system, defined as the mean kinetic 
energy per star divided by the global mean mass (in order to have a "temperature" expressed in square 
velocity units). 

In this plot we show the core collapse situation corresponding to Figs. (6.14) and (6.15): We consider 
a 10 component cluster with the characteristics explained before. The mean temperature is defined like 
L"i ' 7i/£«„ where n\ is the numerical local density for component i. This corresponds to the mean 
kinetic energy per star. We can see in the plot that it is about the same as the heaviest component in 
the inner regions, even though one could think that segregation should not have set in in the beginning. 
This is due to the fact that the moment does not correspond to exactly the initial moment, T = 0. On the 
right, we can already see how the mean central temperature moves back (solid black line) and the most 
massive component (dashed red line) increases. For later times, the kinetic energies of the different 
components rise at the inner part of the cluster and the most massive one approaches the sum of all of 
them. This is even more evident in the last plot, where all components' temperatures have sunken except 
for the corresponding to the most massive one. At the exterior zones, the mean temperature is much 
lower than the most heavy component because the system did not have evolution towards equipartition 
and, so, the velocity dispersion a v of the heavy components is the same that the other's component 
velocity dispersion but their kinetic energy is much larger (in terms of mass). 

At the current moment we are tackling with these and another configurations, but we have to first resolve 
the problem of BH growth and grid resolution, because as the BH grows, the grid gets too close to 
the centre. This problem is being now studied to carry on with new calculations from the cluster re- 
expansion onwards caused by the BH growth. 
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Fig. 6.14: Lagrange radii evolution for a 10 components calculation with a seed BH and stellar remnants 
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Fig. 6.15: Density profiles in a multi-mass system with seed BH before and after core-collapse (see text) 
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Fig. 6.16: Different moments in the evolution of the cluster's temperature for a 10 stellar components 
system with a seed BH 
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eidelberg is a little city in the south of Germany. It has a nice old bridge, a romantic 
castle, a very concurred old town and -I was told- a couple of native inhabitants. I 
spent in this city about five years of my life, and I daresay they were the most intesive 
ones. A lot of things happened. "For instance", I got to know Sabine after one and a 
half weeks of my arrival. Very soon we started the relationship that would conduct us 
in a inavoidably retoric way to marriage. 
If somebody asked me to tell the first thing that comes to my mind when I think about these years, I 
would answer the following... 

sabine bridge potatoes kitchen Schwalbe Jorge Mayer Arger fill de sa mare Francine Vienna bridge cold cold 
cold cold mountains sauerkraut Frank Wixenpa'aslocas wine Oberseminar hin und her fahren Wien cold cold 
cold HauptstraGe Emil potato Mike Deiters 4| [|n4 g3nd m3$ r4r374 Teeminar teaminar teminable @Weberstra6e 
13 Sant Petersburg Masha Sergei Russian Russian Mahos Tilemahos Vienna Sabine Gerhabt-Hauptmann StraBe 
Mutter Umzug moving moving CRACK oh my goodness that was my back Basel Alfonso shop buy smoke Dante 
Inferno work work woman in pink tracksuit yellow umbrella Nickel Nickel what for und viele iuius for a long 
while escape go back to Spain and walk walk walk in the night DanzigerstraBe John Lucie tinto Sabine is back 
al leluia per que peta borinot turbi (A) anxious -xit anxiety dean- and uprooting syndrome quatsch better cold 
dark night Vienna Russian Sprachlabor Marc Russian Silke Eva Marion work work Marc wark Mork Russian 
workian Freitag BISdsinn in der Kiiche paper light see end pictures marriage vacherin gluant Marc's supermongo 
subroutines again Dante Inferno but smoothed end 

This section cannot be short at all because I am in hock to a lot of people. I cannot help feeling very 
grateful to them and they deserve it. On the other hand, being a social animal, as only a genetically 
galician but in Valencia born person can be, it would be a crass error and an unpardonable insult not to 
talk about the friends I made during these years. 

First of all, I would like to thank here my advisor, Rainer Spurzem, for giving me the opportunity of 
coming to Germany to do a Ph.D. He encouraged me always when things were not working as they 
should (a bad habit of the process of becoming a doctor, according to a lot of people). He and Andreas 
Just treated me as if I had been in the group for years from the very beginning. 

The French-speaking commonwealth, Marc Freitag and Francine Leeuwin, carried out a very impor- 
tant role during these years... I daresay Marc is responsible for me becoming a scientist, if you happen 
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to have any complaint, please refer to him. They gave me also something much more valuable than 
their help, namely their friendship. I am very grateful to both of them for a lot of things... Francine, for 
instance, introduced me into the wonderful world of french spirits (I am not talking about the agent or 
subject of vital and spiritual functions, whether spiritual or material). It helped a lot to survive the ab- 
solutely sun-less winter evenings of Heidelberg after a long full working day. I got to know Marc when 
he visited our group in Heidelberg in March 2001. 1 felt as if we had been friends since always from the 
very first beer at the Markt Stiibel (a place which we haunted later). When he moved to Heidelberg we 
became inseparable and even neighbours. Among innumerable things, Marc made me rediscover Edika 
in his full splendour. It would not be an exaggeration at all to say that both of them, Marc and Francine, 
have become a part of my family (in its full southern meaning). 

Elaborating on nice people, I have to mention here two very nice aussies, Bruce and Sheila... I mean 
John and Lucie Benz (nee Chapman), who made me feel "at home" from the very beginning. Every- 
thing began with a red wine, carried on with tinto and jamon (pata negra, of course) and finished with 
two marriages, theirs (on next Saturday, by the way) and mine. 

As regards the Kersche community, it is a must to say how nice it was for me to share a place in the 
Kirchheim community with our fellow-beings Eva, Marion, Silke and -inevitably- Marc. The iberians 
would not be happy if I forgot them here and so I have to address a few words to my two favourite 
tableland people, Jorge Peflarrubia, el Doctor Xungalin, companero infatigable de diversas peripecias 
acontecidas en el consuetudinario devenir, por no llamarlo divagar, del instituto (lease cocina, broncas 
diversas con el aparato administrativo ariense, puerta corrediza y un sinfin mas de cosas) and el lobo 
mesetario David Lopez and parienta, Elise Schieck, for his spirituous southern chocolate, republican 
ideas and Swedish food. And here I cannot -I must not- forget Pablo Gonzalez, Polito, den entspan- 
ntesten Canario der Welt, for his natural sedateness, papas arrugadas con mojo rojo (y verde) and ron 
canario. 

Respect to the ARI people, Gemot Burkhardt's predisposition to make any kind of favour I asked for 
nicely surprised me everytime. He supported us, Jorge and me, in our little "domestic" conflicts at the 
Institute and was always very friendly to us. Somebody like he can make you see things in a different 
way when you start as a newbie and have to confront difficulties and people who are not as nice as he is. 
Christoph Eichhorn, Kristin Warnick, Chingis Omarov, Toshio Tsuchiya, Eliani Ardi, Ikbal 
Aryfanto, Jose Fiestas, Emil Khalisi, Michael Fellhauer, Stefan Deiters, Andreas Ernst, Patrick 
Glaschke, Pavel Kroupa and Christian Boily were very cordial and amiable to me during my stay at 
ARI. 

Ernst Pendl, Eanstl, bin i dangba fia dem Vada sei Dochda Sabine und seine Begeisterung fur As- 
tronomic Im Kosmos umanaundstrawanzn... per aspera ad astra. - Heast! geh gibtsn des a, Spania de 
wos Weanerisch kennan! But yes, there are, and I am one. Da hlb i earns einegsat, net! 

Ich weiB Christine ob der nach einer Woche meines Daseins in Deutschland geschlossenen Bekannt- 
schaft zu ihrer Nachfahrin sehr zu schatzen. 

Margaret Mehls, Margarita war immer auBerst liebenswiirdig zu mir in der WeberstraBe 13 (Aussen- 
stelle vom ARI). Es macht einen riesengroBen Unterschied, ob deine "Nachbarn" (eigentlich Mitbe- 
wohner, nachdem ich fast den ganzen Tag hier bin) nett sind oder nicht. In ihrem Falle (und in dem von 
Herrn Mehls, selbstredend) sollte man eigentlich eher iiber Herzlichkeit als Nettigkeit reden. 

Matt Bonner em va ajudar amb l'Angles de la petita introduccio que vaig fer, car se'm feia molt 
envitricollat de traduir F Homer i les meues paraules a la seua vernacla saxona. Ell, habituat com esta a 
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traduir coses al Catala per a Debian, ho feu en un tres i no-res i, el millor, perque si, sense coneixer-me 
de quasi res, nomes del forum de traduccio. Amb persones com Mateu Bonet (com a ell li agrada que 
li diguen), hom pot creure que es pot recuperar les esperances. No vaig gosar de demanar ningu de 
Micro$oft. Soc optimista, pero no tant. 

Going back a bit in history, just over a couple of years, when I was studying at the Universty of Valencia, 
Vicent Martinez Sancho played a decisive role in a lot a different aspects of my education. He not only 
managed to unfailingly wake up in all his students an enthusiastic interest for the subjects he teaches, but 
also won our trust and affection. Besides Physics, he taught us how to tackle with life, which is a much 
more difficult subject. Putting it into his words, "Esser dotat per a les matematiques no te cap merit. Car 
merit ve de meritum, del participi de merere 'mereixer'. Es aixo, una qualitat natural... Com el beneit 
del cabas del meu poble, que fa trues molt enginyosos amb un palet. El que faces amb la teua vida, com 
l'enfoques i quin tipus d'esser huma esdevingues, aixo si que mereix". His concern and attentiveness 
for culture (understood in all its broad meaning) influenced, influences, me until now. 

Cuando, en los ultimos dfas de la redaction de la tesis, me puse a repasar estas paginas, me di cuenta 
de que el nombre de mi madre no aparecfa en ninguna de ellas. ^Me habfa olvidado de ella? De mi 
madre?! No era posible. En ese momento, en que estaba pensando (contecimiento extraordinario), me 
di cuenta de que habfa sido porque para mi es tan evidente todo lo que le debo -que es, eso, todo- y 
lo mucho que la quiero -que es un muchfsimo mucheante mucheado- que, de hecho, el nombre, no es 
importante. Lo que si es importante es todo lo que ella me ha aportado y aporta en la vida. Gracias, 
Luisa. 

Uber Dich brauche ich wenig zu sagen, Sabine, denn es ist nicht nur fur mich, sondern fur alle klar, wie 
sehr ich Dich liebe. Wir haben ein ganzes Leben vor uns, das gerade anfangt. 

3$ gl4 g'4$f g41 3$m3nd4 31 nOm d31($) $4ny0($) Els Etsus: 3rr|gu|, 3rm|[| | 3rb3r3 d'b4n|4 b3r dOd 
3r gOn rOdlll g4'n$ h4m blrd4d 4gu3$d$ 4ny$. 

GOn d | -D-4rl ! ! ! 1 (sorry for that, but it was also a must). 
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BuUUaimf I think this is really the end. I do not know what else I could tell 
you, after 16rpages of wanton hullabaloo... Now it is time to have a carajillo... 
Let me pythonisf it... 



This work was written with BT[=X (what else?) and processed with a Debian Sid- based system, VINDOBONA 



